A critical evaluation of different QM/MM frontier treatments with SCC-DFTB as the QM method

The performance of different link atom based frontier treatments in QM/MM simulations was evaluated critically with SCC-DFTB as the QM method. In addition to the analysis of gas-phase molecules as in previous studies, an important element of the present work is that chemical reactions in realistic enzyme systems were also examined. The schemes tested include all options available in the program CHARMM for SCC-DFTB/MM simulation, which treat electrostatic interactions due to the MM atoms close to the QM/MM boundary in different ways. In addition, a new approach, the divided frontier charge (DIV), has been implemented in which the partial charge associated with the frontier MM atom ("link host") is evenly distributed to the other MM atoms in the same group. The performance of these schemes was evaluated based on properties including proton affinities, deprotonation energies, dipole moments, and energetics of proton transfer reactions. Similar to previous work, it was found that calculated proton affinities and deprotonation energies of alcohols, carbonic acids, amino acids, and model DNA bases are very sensitive to the link atom scheme; the commonly used single link atom approach often gives error on the order of 15 to 20 kcal/mol. Other schemes give better and, on average, mutually comparable results. For proton transfer reactions, encouragingly, both activation barriers and reaction energies are fairly insensitive (within a typical range of 2-4 kcal/mol) to the link atom scheme due to error cancellation, and this was observed for both gas-phase and enzyme systems. Therefore, the effect of using different link atom schemes in QM/MM simulations is rather small for chemical reactions that conserve the total charge. Although the current study used an approximate DFT method as the QM level, the observed trends are expected to be applicable to QM/MM methods with use of other QM approaches. This observation does not mean to encourage QM/MM simulations without careful benchmark in the study of specific systems, rather it emphasizes that other technical details, such as the treatment of long-range electrostatics, tend to play a more important role and need to be handled carefully.     

Importance of van der Waals Interactions in QM/MM Simulations

The importance of accurately treating van der Waals interactions between the quantum mechanical (QM) and molecular mechanical (MM) atoms in hybrid QM/MM simulations has been investigated systematically. First, a set of van der Waals (vdW) parameters was optimized for an approximate density functional method, the self-consistent charge-tight binding density functional (SCC-DFTB) approach, based on small hydrogen-bonding clusters. The sensitivity of condensed phase observables to the SCC-DFTB vdW parameters was then quantitatively investigated by SCC-DFTB/MM simulations of several model systems using the optimized set and two sets of extreme vdW parameters selected from the CHARMM22 forcefield. The model systems include a model FAD molecule in solution and a solvated enediolate, and the properties studied include the radial distribution functions of water molecules around the solute (model FAD and enediolate), the reduction potential of the model FAD and the potential of mean force for an intramolecular proton transfer in the enediolate. Although there are noticeable differences between parameter sets for gas-phase clusters and solvent structures around the solute, thermodynamic quantities in the condensed phase (e.g., reduction potential and potential of mean force) were found to be less sensitive to the numerical values of vdW parameters. The differences between SCC-DFTB/MM results with the three vdW parameter sets for SCC-DFTB atoms were explained in terms of the effects of the parameter set on solvation. The current study has made it clear that efforts in improving the reliability of QM/MM methods for energetical properties in the condensed phase should focus on components other than van der Waals interactions between QM and MM atoms.     

Reliable treatment of electrostatics in combined QM/MM simulation of macromolecules

A robust approach for dealing with electrostatic interactions for spherical boundary conditions has been implemented in the QM/MM framework. The development was based on the generalized solvent boundary potential (GSBP) method proposed by Im et al. [J. Chem. Phys. 114, 2924 (2001)], and the specific implementation was applied to the self-consistent-charge density-functional tight-binding approach as the quantum mechanics (QM) level, although extension to other QM methods is straightforward. Compared to the popular stochastic boundary-condition scheme, the new protocol offers a balanced treatment between quantum mechanics/molecular mechanics (QM/MM) and MM/MM interactions; it also includes the effect of the bulk solvent and macromolecule atoms outside of the microscopic region at the Poisson-Boltzmann level. The new method was illustrated with application to the enzyme human carbonic anhydrase II and compared to stochastic boundary-condition simulations using different electrostatic treatments. The GSBP-based QM/MM simulations were most consistent with available experimental data, while conventional stochastic boundary simulations yielded various artifacts depending on different electrostatic models. The results highlight the importance of carefully treating electrostatics in QM/MM simulations of biomolecules and suggest that the commonly used truncation schemes should be avoided in QM/MM simulations, especially in simulations that involve extensive conformational samplings. The development of the GSBP-based QM/MM protocol has opened up the exciting possibility of studying chemical events in very complex biomolecular systems in a multiscale framework.     

pKa analysis for the zinc-bound water in human carbonic anhydrase II: Benchmark for "multiscale" QM/MM simulations and mechanistic implications

To quantitatively explore the applicability of the generalized solvent boundary potential (GSBP) based QM/MM approach as a "multiscale" framework for studying chemical reactions in biomolecules, the structural and energetic properties of the Human Carbonic Anhydrase II (CAII) are analyzed and compared to those from periodic boundary condition (PBC) simulations and available experimental data. Although the atomic fluctuations from GSBP based simulations are consistently lower compared to those from PBC simulations or crystallographic data, the fluctuations and internal coordinate distributions for residues in the proximity of the active site as well as diffusion constants of active-site water molecules are fairly well described by GSBP simulations. The pKa of the zinc-bound water, calculated with a SCC-DFTB/MM-GSBP based thermodynamic integration approach, agrees well with experiments for the wild type CAII. For the E106Q mutant, however, a 9 pKa unit downward shift relative to the wild type is found in contrast with previous experiments that found little change. This dramatic discrepancy signals a possible change in the mechanism for the interconversion between CO2/HCO3- in the E106Q mutant, which may be similar to the bicarbonate mediated mechanism proposed for the Co2+ substituted CAII (J. Am. Chem. Soc. 2001, 123, 5861).1 The study highlights pKa analyses as a valuable approach for quantitatively validating the computational model for complex biomolecules as well as for revealing energetic properties intimately related to the chemical process of interest.     

The implementation of a fast and accurate QM/MM potential method in Amber

Version 9 of the Amber simulation programs includes a new semi-empirical hybrid QM/MM functionality. This includes support for implicit solvent (generalized Born) and for periodic explicit solvent simulations using a newly developed QM/MM implementation of the particle mesh Ewald (PME) method. The code provides sufficiently accurate gradients to run constant energy QM/MM MD simulations for many nanoseconds. The link atom approach used for treating the QM/MM boundary shows improved performance, and the user interface has been rewritten to bring the format into line with classical MD simulations. Support is provided for the PM3, PDDG/PM3, PM3CARB1, AM1, MNDO, and PDDG/MNDO semi-empirical Hamiltonians as well as the self-consistent charge density functional tight binding (SCC-DFTB) method. Performance has been improved to the point where using QM/MM, for a QM system of 71 atoms within an explicitly solvated protein using periodic boundaries and PME requires less than twice the cpu time of the corresponding classical simulation.     

51V NMR chemical shifts calculated from QM/MM models of vanadium chloroperoxidase

(51)V NMR chemical shifts calculated from QM/MM-optimized (QM=quantum mechanical; MM=molecular mechanical) models of vanadium-dependent chloroperoxidase (VCPO) are presented. An extensive number of protonation states for the vanadium cofactor (active site of the protein) and a number of probable positional isomers for each of the protonation states are considered. The size of the QM region is increased incrementally to observe the convergence behavior of the (51)V NMR chemical shifts. A total of 40 models are assessed by comparison to experimental solid-state (51)V NMR results recently reported in the literature. Isotropic chemical shifts are found to be a poor indicator of the protonation state; however, anisotropic chemical shifts and the nuclear quadrupole tensors appear to be sensitive to changes in the proton environment of the vanadium nuclei. This detailed investigation of the (51)V NMR chemical shifts computed from QM/MM models provides further evidence that the ground state is either a triply protonated (one axial water and one equatorial hydroxyl group) or a doubly protonated vanadate moiety in VCPO. Particular attention is given to the electrostatic and geometric effects of the protein environment. This is the first study to compute anisotropic NMR chemical shifts from QM/MM models of an active metalloprotein for direct comparison with solid-state MAS NMR data. This theoretical approach enhances the potential use of experimental solid-state NMR spectroscopy for the structural determination of metalloproteins.     

MM and QM/MM Modeling of Threonyl-tRNA Synthetase: Model Testing and Simulations

Aminoacyl-tRNA synthetases are centrally important enzymes in protein synthesis. We have investigated threonyl-tRNA synthetase from E. coli, complexed with reactants, using molecular mechanics and combined quantum mechanical/molecular mechanical (QM/MM) techniques. These modeling methods have the potential to provide molecular level understanding of enzyme catalytic processes. Modeling of this enzyme presents a number of challenges. The procedure of system preparation and testing is described in detail. For example, the number of metal ions at the active site, and their positions, were investigated. Molecular dynamics simulations suggest that the system is most stable when it contains only one magnesium ion, and the zinc ion is removed. Two different QM/MM methods were tested in models based on the findings of MM molecular dynamics simulations. AM1/CHARMM calculations resulted in unrealistic structures for the phosphates in this system. This is apparently due to an error of AM1. PM3/CHARMM calculations proved to be more suitable for this enzyme system. These results will provide a useful basis for future modeling investigations of the enzyme mechanism and dynamics.     

What is so special about Arg 55 in the catalysis of cyclophilin A? insights from hybrid QM/MM simulations

Potential of mean force (PMF) simulations with a hybrid QM/MM potential function were used to analyze the catalytic mechanism of human cyclophilin A (CypA). PMF calculations were performed for proline isomerization of peptides in solution, the wild-type CypA, and several CypA mutants. With an approximate density functional theory, the self-consistent-charge density functional tight binding (SCC-DFTB) as the QM level, and CHARMM 22 force field as MM, satisfactory energetics compared to available experiments were obtained. Calculations for the Arg55Ala and zero-charge-Arg55 mutants clearly indicated that Arg 55 significantly stabilizes the isomerization transition state through electrostatic interactions. However, the decrease in the average distance (thus the increase in interaction) between Arg 55 and the substrate amide N in going from the stable states to the transition state is mainly due to the pyramidalization of the amide N rather than motions associated with Arg 55. Although the nanosecond simulations cannot exclude the existence of sub-millisecond collective motions proposed on the basis of recent elegant NMR relaxation and line-shape analyses, the energetics obtained for the various enzyme systems here indicate that the contribution from motions of active site residues to catalysis is expected to be small. Instead, the present simulations support that the structural stability rather than mobility of the preorganized active site is more important. Through hydrogen-bonding interactions among the substrate, Arg 55, Gln 63, and Asn 102, the active site of the wild-type enzyme is structurally very stable and puts Arg 55 in a favorable position to perform its catalytic role in the transition state. This is further illustrated with the somewhat unexpected prediction that Arg55Lys is largely catalytically inactive, because Lys does not have the unique bifurcating construct of the guanidino group in Arg and thus the active site of Arg55Lys cannot accommodate Lys in a position capable of providing electrostatic stabilization of the isomerization transition state. Among all the enzyme systems studied, the wild-type CypA is the only one that selects the syn/exo transition state, while the syn/endo conformation is also present in the mutants, which is another reason for their higher barriers. Finally, the present analysis indicated that the population of near-attack-conformations (NAC) is not relevant to catalysis in CypA.     

Interactions between phosphate and water in solution: a natural bond orbital based analysis in a QM/MM framework

The natural bond orbital (NBO) and natural energy decomposition analysis (NEDA) calculations are used to analyze the interaction between mono-methyl phosphate-ester (MMP) and its solvation environment in a combined quantum mechanical/molecular mechanical (QM/MM) framework. The solute-solvent configurations are generated using a specific parametrization of the self-consistent-charge density functional tight-binding (SCC-DFTB) model for the MMP and TIP3P for water. The NBO and NEDA calculations are done with several QM/MM partitioning schemes with HF/6-31+G** as the QM level. Regardless of the size of the QM region, a notable amount of charge transfer is observed between MMP and the neighboring water molecules and the charge-transfer interactions are, in the NEDA framework, as important as the electric (electrostatic and polarization) components. This work illustrates that NBO based analyses are effective tools for probing intermolecular interactions in condensed phase systems.     

Mechanism of C-terminal intein cleavage in protein splicing from QM/MM molecular dynamics simulations

Protein splicing is a post-translational process in which a biologically inactive protein is activated by the release of a segment denoted as an intein. The process involves four steps. In the third, the scission of the intein takes place after the cyclization of the last amino acid of the segment, an asparagine. Little is known about the chemical reaction necessary for this cyclization. Experiments demonstrate that two histidines (the penultimate amino acid of the intein, and a histidine located 10 amino acids upstream) are relevant in the cyclization of the asparagine. We have investigated the mechanism and determinants of reaction in the GyrA intein focusing on the requirements for asparagine activation for its cyclization. First, the influence that the protonation states of these two histidines have on the orientation of the asparagine side chain is investigated by means of molecular dynamics simulation. Molecular dynamics simulations using the CHARMM27 force field were carried out on the three possible protonation states for each of these two histidines. The results indicate that the only protonation state in which the conformation of the system is suitable for cyclization is when the penultimate histidine is fully protonated (positively charged), and the upstream histidine is in the His(ε) neutral tautomeric form. The free energy profile for the reaction in which the asparagine is activated by a proton transfer to the upstream histidine is presented, computed by hybrid quantum mechanics/molecular mechanics (QM/MM) umbrella sampling molecular dynamics at the SCCDFTB/CHARMM27 level of theory. The calculated free energy barrier for the reaction is 19.0 kcal mol(-1). B3LYP/6-31+G(d) QM/MM single-point calculations give a qualitatively a similar energy profile, although with somewhat higher energy barriers, in good agreement with the value derived from experiment of 25 kcal mol(-1) at 60 °C. QM/MM molecular dynamics simulations of the reactant, activated reactant and intermediate states highlight the importance of the Arg181-Val182-Asp183 segment in catalysing the reaction. Overall, the results indicate that nucleophilic activation of the asparagine for its cyclization by the upstream histidine acting as the base is a plausible mechanism for the C-terminal cleavage in protein splicing.     

How does the cAMP-dependent protein kinase catalyze the phosphorylation reaction: an ab initio QM/MM study

We have carried out density functional theory QM/MM calculations on the catalytic subunit of cAMP-dependent protein kinase (PKA). The QM/MM calculations indicate that the phosphorylation reaction catalyzed by PKA is mainly dissociative, and Asp166 serves as the catalytic base to accept the proton delivered by the substrate peptide. Among the key interactions in the active site, the Mg(2+) ions, glycine rich loop, and Lys72 are found to stabilize the transition state through electrostatic interactions. On the other hand, Lys168, Asn171, Asp184, and the conserved waters bound to Mg(2+) ions do not directly contribute to lower the energy barrier of the phosphorylation reaction, and possible roles for these residues are proposed. The QM/MM calculations with different QM/MM partition schemes or different initial structures yield consistent results. In addition, we have carried out 12 ns molecular dynamics simulations on both wild type and K168A mutated PKA, respectively, to demonstrate that the catalytic role of Lys168 is to keep ATP and substrate peptide in the near-attack reactive conformation.     

Exploring SCC-DFTB paths for mapping QM/MM reaction mechanisms

A new first-order procedure for locating transition structures (TS) that employs hybrid quantum mechanical/molecular mechanical (QM/MM) potentials has been developed. This new technique (RPATh+RESD) combines the replica path method (RPATh) and standard reaction coordinate driving (RCD) techniques in an approach that both efficiently determines reaction barriers and successfully eliminates two key weaknesses of RCD calculations (i.e., hysteresis/discontinuities in the path and the sequential nature of the RCD procedure). In addition, we have extended CHARMM's QM/MM reaction pathway methods, the RPATh and nudged elastic band (NEB) methods, to incorporate SCC-DFTB wave functions. This newly added functionality has been applied to the chorismate mutase-catalyzed interconversion of chorismate to prephenate, which is a key step in the shikimate pathway of bacteria, fungi, and other higher plants. The RPATh+RESD barrier height (DeltaE=5.7 kcal/mol) is in good agreement with previous results from full-energy surface mapping studies (Zhang, X.; Zhang, X.; Bruice, T. C. Biochemistry 2005, 44, 10443-10448). Full reaction paths were independently mapped with RPATh and NEB methods and showed good agreement with the final transition state from the RPATh+RESD "gold standard" and previous high-level QM/MM transition states (Woodcock, H. L.; Hodoscek, M.; Gilbert, T. B.; Gill, P. M. W.; Schaefer, H. F.; Brooks, B. R. J. Comput. Chem. 2007, 28, 1485-1502). The SCC-DFTB TS geometry most closely approximates the MP2/6-31+G(d) QM/MM result. However, the barrier height is underestimated and possibly points to an area for improvement in SCC-DFTB parametrization. In addition, the steepest descents (SD) minimizer for the NEB method was modified to uncouple the in-path and off-path degrees of freedom during the minimization, which significantly improved performance. The convergence behavior of the RPATh and NEB was examined for SCC-DFTB wave functions, and it was determined that, in general, both methods converge at about the same rate, although the techniques used for convergence may be different. For instance, RPATh can effectively use the adopted basis Newton-Raphson (ABNR) minimizer, where NEB seems to require a combination of SD and ABNR.     

New QM/MM implementation of the DFTB3 method in the gromacs package

The approximate density‐functional tight‐binding theory method DFTB3 has been implemented in the quantum mechanics/molecular mechanics (QM/MM) framework of the Gromacs molecular simulation package. We show that the efficient smooth particle–mesh Ewald implementation of Gromacs extends to the calculation of QM/MM electrostatic interactions. Further, we make use of the various free‐energy functionalities provided by Gromacs and the PLUMED plugin. We exploit the versatility and performance of the current framework in three typical applications of QM/MM methods to solve biophysical problems: (i) ultrafast proton transfer in malonaldehyde, (ii) conformation of the alanine dipeptide, and (iii) electron‐induced repair of a DNA lesion. Also discussed is the further development of the framework, regarding mostly the options for parallelization. © 2015 Wiley Periodicals, Inc.     

Fragmentation-based QM/MM simulations: length dependence of chain dynamics and hydrogen bonding of polyethylene oxide and polyethylene in aqueous solutions

Molecular fragmentation quantum mechanics (QM) calculations have been combined with molecular mechanics (MM) to construct the fragmentation QM/MM method for simulations of dilute solutions of macromolecules. We adopt the electrostatics embedding QM/MM model, where the low-cost generalized energy-based fragmentation calculations are employed for the QM part. Conformation energy calculations, geometry optimizations, and Born-Oppenheimer molecular dynamics simulations of poly(ethylene oxide), PEO(n) (n = 6-20), and polyethylene, PE(n) ( n = 9-30), in aqueous solution have been performed within the framework of both fragmentation and conventional QM/MM methods. The intermolecular hydrogen bonding and chain configurations obtained from the fragmentation QM/MM simulations are consistent with the conventional QM/MM method. The length dependence of chain conformations and dynamics of PEO and PE oligomers in aqueous solutions is also investigated through the fragmentation QM/MM molecular dynamics simulations.     

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.     

Spectral tuning in visual pigments: an ONIOM(QM:MM) study on bovine rhodopsin and its mutants

We have investigated geometries and excitation energies of bovine rhodopsin and some of its mutants by hybrid quantum mechanical/molecular mechanical (QM/MM) calculations in ONIOM scheme, employing B3LYP and BLYP density functionals as well as DFTB method for the QM part and AMBER force field for the MM part. QM/MM geometries of the protonated Schiff-base 11- cis-retinal with B3LYP and DFTB are very similar to each other. TD-B3LYP/MM excitation energy calculations reproduce the experimental absorption maximum of 500 nm in the presence of native rhodopsin environment and predict spectral shifts due to mutations within 10 nm, whereas TD-BLYP/MM excitation energies have red-shift error of at least 50 nm. In the wild-type rhodopsin, Glu113 shifts the first excitation energy to blue and accounts for most of the shift found. Other amino acids individually contribute to the first excitation energy but their net effect is small. The electronic polarization effect is essential for reproducing experimental bond length alternation along the polyene chain in protonated Schiff-base retinal, which correlates with the computed first excitation energy. It also corrects the excitation energies and spectral shifts in mutants, more effectively for deprotonated Schiff-base retinal than for the protonated form. The protonation state and conformation of mutated residues affect electronic spectrum significantly. The present QM/MM calculations estimate not only the experimental excitation energies but also the source of spectral shifts in mutants.     

Combined QM/MM study of the opsin shift in bacteriorhodopsin

Combined quantum mechanical and molecular mechanical (QM/MM) calculations and molecular dynamics simulations of bacteriorhodopsin (bR) in the membrane matrix have been carried out to determine the factors that make significant contributions to the opsin shift. We found that both solvation and interactions with the protein significantly shifts the absorption maximum of the retinal protonated Schiff base, but the effects are much more pronounced in polar solvents such as methanol, acetonitrile, and water than in the protein environment. The differential solvatochromic shifts of PSB in methanol and in bR leads to a bathochromic shift of about 1800 cm(-1). Because the combined QM/MM configuration interaction calculation is essentially a point charge model, this contribution is attributed to the extended point-charge model of Honig and Nakanishi. The incorporation of retinal in bR is accompanied by a change in retinal conformation from the 6-s-cis form in solution to the 6-s-trans configuration in bR. The extension of the pi-conjugated system further increases the red-shift by 2400 cm(-1). The remaining factors are due to the change in dispersion interactions. Using an estimate of about 1000 cm(-1) in the dispersion contribution by Houjou et al., we obtained a theoretical opsin shift of 5200 cm(-1) in bR, which is in excellent agreement with the experimental value of 5100 cm(-1). Structural analysis of the PSB binding site revealed the specific interactions that make contributions to the observed opsin shift. The combined QM/MM method used in the present study provides an opportunity to accurately model the photoisomerization and proton transfer reactions in bR.     

The CHARMM–TURBOMOLE interface for efficient and accurate QM/MM molecular dynamics,free energies,and excited state properties

The quantum mechanical (QM)/molecular mechanical (MM) interface between Chemistry at HARvard Molecular Mechanics (CHARMM) and TURBOMOLE is described. CHARMM provides an extensive set of simulation algorithms, like molecular dynamics (MD) and free energy perturbation, and support for mature nonpolarizable and Drude polarizable force fields. TURBOMOLE provides fast QM calculations using density functional theory or wave function methods and excited state properties. CHARMM–TURBOMOLE is well‐suited for extended QM/MM MD simulations using first principles methods with large (triple‐ζ) basis sets. We demonstrate these capabilities with a QM/MM simulation of Mg²⁺(aq), where the MM outer sphere water molecules are represented using the SWM4‐NDP Drude polarizable force field and the ion and inner coordination sphere are represented using QM PBE, PBE0, and MP2 methods. The relative solvation free energies of Mg²⁺ and Zn²⁺ were calculated using thermodynamic integration. We also demonstrate the features for excited state properties. We calculate the time‐averaged solution absorption spectrum of indole, the emission spectrum of the indole excited state, and the electronic circular dichroism spectrum of an oxacepham. © 2014 Wiley Periodicals, Inc.     

Interfacing Q-Chem and CHARMM to perform QM/MM reaction path calculations

A hybrid quantum mechanical/molecular mechanical (QM/MM) potential energy function with Hartree-Fock, density functional theory (DFT), and post-HF (RIMP2, MP2, CCSD) capability has been implemented in the CHARMM and Q-Chem software packages. In addition, we have modified CHARMM and Q-Chem to take advantage of the newly introduced replica path and the nudged elastic band methods, which are powerful techniques for studying reaction pathways in a highly parallel (i.e., parallel/parallel) fashion, with each pathway point being distributed to a different node of a large cluster. To test our implementation, a series of systems were studied and comparisons were made to both full QM calculations and previous QM/MM studies and experiments. For instance, the differences between HF, DFT, MP2, and CCSD QM/MM calculations of H2O...H2O, H2O...Na+, and H2O...Cl- complexes have been explored. Furthermore, the recently implemented polarizable Drude water model was used to make comparisons to the popular TIP3P and TIP4P water models for doing QM/MM calculations. We have also computed the energetic profile of the chorismate mutase catalyzed Claisen rearrangement at various QM/MM levels of theory and have compared the results with previous studies. Our best estimate for the activation energy is 8.20 kcal/mol and for the reaction energy is -23.1 kcal/mol, both calculated at the MP2/6-31+G(d)//MP2/6-31+G(d)/C22 level of theory.