MM simulations of Henry and methyl transfer reactions.

HHS Public Access Author manuscript Author Manuscript

J Comput Chem. Author manuscript; available in PMC 2016 October 15. Published in final edited form as: J Comput Chem. 2015 October 15; 36(27): 2064–2074. doi:10.1002/jcc.24045.

Application of a BOSS – Gaussian Interface for QM/MM Simulations of Henry and Methyl Transfer Reactions Jonah Z. Vilseck, Jakub Kostal, Julian Tirado-Rives, and William L. Jorgensen* Department of Chemistry, Yale University, New Haven, CT 06520-8107USA

Abstract Author Manuscript Author Manuscript

Hybrid quantum mechanics and molecular mechanics (QM/MM) computer simulations have become an indispensable tool for studying chemical and biological phenomena for systems too large to treat with quantum mechanics alone. For several decades, semi-empirical QM methods have been used in QM/MM simulations. However, with increased computational resources, the introduction of ab initio and density function methods into on-the-fly QM/MM simulations is being increasingly preferred. This adaptation can be accomplished with a program interface that tethers independent QM and MM software packages. This report introduces such an interface for the BOSS and Gaussian programs, featuring modification of BOSS to request QM energies and partial atomic charges from Gaussian. A customizable C-shell linker script facilitates the interprogram communication. The BOSS–Gaussian interface also provides convenient access to Charge Model 5 (CM5) partial atomic charges for multiple purposes including QM/MM studies of reactions. In this report, the BOSS–Gaussian interface is applied to a nitroaldol (Henry) reaction and two methyl transfer reactions in aqueous solution. Improved agreement with experiment is found by determining free-energy surfaces with MP2/CM5 QM/MM simulations than previously reported investigations employing semiempirical methods.

Keywords QM/MM; BOSS; Gaussian 09; Charge Model 5; Henry reaction; methyl transfer reaction; reaction path

Introduction Author Manuscript

From its introduction in the 1970s[1] to its recognition by the 2013 Nobel Prize in Chemistry,[2] hybrid quantum and molecular mechanics (QM/MM) computer models have proven indispensable for describing molecular interactions in every field of chemistry. Organic reactions, enzyme catalysis, solvation, and protein–ligand binding represent only a few popular applications.[3,4] The appeal is clear. QM/MM calculations combine the ability of modern day force fields and statistical mechanics to model systems containing hundreds

Correspondence to: William L. Jorgensen ([email protected]). Author Contributions: All authors conceived and designed experiments. JZV and JK programmed the interface, and performed the QM/MM simulations. All authors analyzed simulation results. JTR and WLJ provided mentoring and support. All authors reviewed and contributed to the manuscript. The authors declare no competing financial interests. **Additional Supporting Information may be found in the online version of this article.

Vilseck et al.

Page 2

Author Manuscript

of thousands of atoms with quantum mechanics’ ability to accurately describe inherent interactions and bond–breaking and bond–forming events. This balance is achieved by dividing a system into two regions. A central or core site of chemical activity is described by QM and all surrounding solute and solvent molecules are described classically with MM. In this manner, typical QM and MM limitations are circumvented, allowing large chemical or biological systems that would be too expensive to model with quantum mechanics alone to be affordably studied.[1–4]

Author Manuscript

Integrating quantum mechanics calculations into ongoing Monte Carlo (MC) or molecular dynamics (MD) statistical mechanics simulations, often referred to as on-the-fly QM/MM, has become a standard practice for describing solution and enzymatic reactions.[3,4] This method allows the QM region to respond and fluctuate with its environment, a useful characteristic for investigating solvent effects on simple organic reactions.[4] To obtain sufficiently converged thermodynamic ensemble averages, however, extensive sampling is required, often necessitating millions or billions of QM calculations.[5] Due to their speed and ease of implementation with various MM software packages, semi-empirical quantum mechanical (SQM) methods, such as AM1,[6] PM3,[7] or PDDG-PM3,[8] have been frequently used.[4] Nevertheless, the limitations of such methods are well documented.[9] In short, SQM methods are normally parameterized for a limited range of elements to reproduce ground state properties and may not describe well systems too different from their training sets, such as transition states.[10–13]

Author Manuscript

With recent growth in computing resources, more accurate ab initio or density functional methods are being increasingly favored for use in on-the-fly QM/MM simulations.[10–20] This can be efficiently accomplished via a program interface, whereby MM and QM programs are coupled together. To date, several interfaces have been reported.[10–20] For example, combinations of AMBER,[21] GROMOS,[22] CHARMM,[23] and TINKER[24] molecular mechanics packages with GAMESS,[25] Gaussian,[26] TURBOMOLE,[27] QChem,[28] and TeraChem[29] quantum mechanics packages have been reported. In general, communication between programs occurs via system calls and file sharing, whereby data is written to disk by one program and read by another program. QM/MM routines already found within some MM packages can be modified to facilitate an interface with ab initio or density functional containing programs. Thus, an MM package is able to retain its original functionality, while extending its ability to request QM energies, forces, and atomic charges, as needed.

Author Manuscript

In this report, an interface between the BOSS and Gaussian programs is described. Gaussian 09 is a popular electronic structure program featuring a wide variety of ab initio and density functional methods.[26] BOSS is a general purpose molecular mechanics engine capable of performing various calculations, such as energy minimizations, normal mode analysis, and conformational searching, in addition to Metropolis Monte Carlo statistical mechanics simulations in both canonical (NVT) and isothermal-isobaric (NPT) ensembles.[5] The OPLS-AA force field is used for classical energy evaluations,[5,30,31] and AM1, PM3, and PDDG/PM3 SQM methods are available for QM and QM/MM calculations. Recent MC applications with BOSS include computing solvation free energies,[32,33] partial molar

J Comput Chem. Author manuscript; available in PMC 2016 October 15.

Vilseck et al.

Page 3

Author Manuscript

volumes of organic solutes,[34] investigation of halogen-bonding effects,[35–37] and mechanistic investigations of reactions in solution.[38–42]

Author Manuscript

Development of a BOSS–Gaussian interface not only provides convenient access to ab initio or density functional methods, but also initiates the ability to incorporate Charge Model 5 (CM5) charges[43] into QM/MM studies or other non-QM/MM MC simulations performed with BOSS. CM5 charges are an attractive alternative to the CM1A or CM3P charge models currently available in BOSS.[44–46] CM5 charges have been parameterized for all atoms in the periodic table and are consistent with a single parameter set over a large range of electronic structure methods and basis sets. This greatly improves the applicability of the model and facilitates its use in a variety of QM calculations. CM5 charges show less conformational dependence and better describe buried atoms in large molecules. Available for use with neutral or ionic molecules, computed CM5 dipole moments are more accurate than those from Hirshfeld, Löwdin, or Mulliken population analyses and are comparable to dipole moments derived from electrostatic potential (ESP) calculations.[43] Recently, CM5 partial atomic charges have been shown to perform well in OPLS-AA based condensed phase simulations.[33] When scaled by a factor of 1.20 for neutral molecules, CM5 partial atomic charges yield mean unsigned errors of ca. 1.0 kcal/mol for computed free energies of hydration and pure liquid heats of vaporization.[33] CM5 charges are also expected to be valuable for use in hybrid QM/MM reaction studies.

Author Manuscript

The BOSS–Gaussian interface discussed here builds upon the development of an earlier interface with GAMESS.[10] The present BOSS–Gaussian interface incorporates key findings from that work, the versatility of CM5 charges, and improved user customizability. This report begins by describing the implementation of the BOSS–Gaussian interface in the context of running QM/MM simulations. The interface is then applied to study three reactions is aqueous solution: a Henry reaction, and two methyl transfer reactions.

The BOSS–Gaussian Interface QM/MM in BOSS The QM/MM routine implemented in BOSS is briefly described.[4,5,47] The chemical system is divided into two regions. Reacting solute(s) or residues of macromolecules are described quantum mechanically. All surrounding explicit solvent molecules, non-reacting solute molecules, or spectator residues are described classically with the OPLS-AA force field.[30,31] Explicit, non-aqueous solvents are represented with the OPLS-AA force field and the TIP4P model is used for water.[48] Metropolis MC statistical mechanics is used for statistical sampling of the entire chemical system (QM/MM/MC).[49]

Author Manuscript

(1)

The total energy of the system is expressed as a sum of QM, MM, and QM/MM energies (eq. 1).[4,5,10,47] In the simplest scenario, no covalent bonds bridge the QM – MM regions and EQM/MM is composed of only nonbonded electrostatic and van der Waals interactions between QM and MM molecules (eq. 2); qi and qj represent atom centered QM and MM


Vilseck et al.

Page 4

Author Manuscript

partial atomic charges, respectively, while σij and εij are standard OPLS-AA Lennard-Jones parameters. Alternatively, if covalent bonds bridge QM – MM regions, such as in an enzyme, the link atom approach can be used.[50] In this instance, EQM/MM also includes appropriate bond stretching, angle bending, and torsional terms for atom pairs crossing the QM – MM boundary.

(2)

Author Manuscript

At present, mechanical embedding is supported for QM/MM simulations in BOSS.[47] Thus, EQM and QM partial atomic charges are derived in the absence of an electrostatic field of surrounding MM point charges, i.e. in vacuum. This is in contrast to electrostatic embedding, which allows QM charges to be explicitly polarized by including MM point charges in the QM calculation. Instead, solute polarization is accounted for in an average sense by scaling the QM charges. Scale factors of 1.1–1.2 have generally been found to reproduce well free energies of hydration of organic molecules, with mean unsigned errors of ca. 1.0 kcal/mol.[32,33,46,47,51] The preferred method to generate QM charges in BOSS has been to use the Cramer and Truhlar Charge Models (CMx), such as CM1A, CM3P, and CM5.[43–45] These methods are optimized to reproduce gas phase dipole moments with errors less than 0.30 D, and they have shown excellent performance in OPLS-AA based calculations.[33,46] It should be noted that charge scaling is only performed for neutral solutes.

Author Manuscript

Although Gaussian 09 supports both mechanical and electrostatic embedding schemes, only mechanical embedding was pursued for interface development at present. In general, mechanical embedding is simpler, less computationally expensive, and presents the most direct route to employ CMx charges.[3–5,38–42,47] Direct polarization of QM partial atomic charges is an appealing advantage of the electrostatic embedding scheme; however, CMx charges would only be affected indirectly since these charges are derived by empirically reweighting Hirshfeld,[43,56] Löwdin,[45] or Mulliken[44] charges. Previous tests to implicitly polarize CM5 charges with the CPCM method indicated that charges remained underpolarized for use in condensed phase simulations. Thus charge scaling would still be recommended for neutral solutes,[33] negating much of the advantage offered by the electrostatic embedding scheme. Furthermore, because QM calculations are not performed for every MC configuration, inclusion of all solute-solvent electrostatic interactions in the QM calculation is not practical.

Author Manuscript

The BOSS–Gaussian Interface Two general strategies have been pursed in previous interface developments.[10–20] Either an MM program has been expanded to incorporate QM calculations or external scripts or programs have been developed to mediate information passing between QM and MM programs. The later strategy was chosen to develop the BOSS–Gaussian interface. BOSS was minimally modified; Gaussian 09 remains unchanged. A C shell linker script, “gausslink”, was written to facilitate communication between the QM and MM programs.[10] Use of a linker script provides flexibility to the interface and reduces


Vilseck et al.

Page 5

Author Manuscript

necessary MM source code changes to a minimum. Thus by simply updating the linker script the interface can be readily modified to function with a variety of QM software packages or updates.[10] To increase the transferability and customizability of the interface, all QM input specifications, such as Gaussian’s %link and %route keyword designations, are located within the linker script itself. Common entries include the name of the checkpoint file, the number of processors to use in parallel, the amount of memory requested, the level of chemical theory to use, the population analysis to use, the title of the calculation, and system charge and multiplicity.

Author Manuscript

Scheme 1 describes the cyclic workflow of the interface. Operations performed by each program are indicated by color: BOSS (blue), gausslink (green), and Gaussian (orange). Upon MC sampling of the QM region, atomic numbers and Cartesian coordinates of all atoms within the QM region are written to disk and gausslink in called. QM coordinates are read by the linker script and the Gaussian input file is assembled. Gaussian is then executed and performs a single-point calculation and a population analysis. Upon termination, the Gaussian output file is parsed by gausslink. EQM and all QM partial atomic charges are written to disk and accessed by BOSS. In BOSS, Etotal is calculated, the Metropolis acceptance test is performed, and MC sampling continues until another MC move of the QM solute is attempted. At that point, the interface workflow repeats. In comparison to noninterface calculations, steps 1–6 in Scheme 1 are replaced by a single semiempirical QM single-point calculation and CMx charge determination within BOSS (SQM arrows in Scheme 1).

Author Manuscript

At present, the following electronic structure methods are scripted into the interface: HF,[52] B3LYP,[53] M06-2X,[54] and MP2.[55] Any basis set available in Gaussian in conjunction with these methods is allowed. The following population analysis methods are also recognized: CM5,[43] Hirshfeld,[56] MK (ESP),[57] and CHelpG.[58] Collectively these lists are limited in comparison to the many options available within Gaussian, but represent a balance of accuracy and speed necessary when performing millions of QM single-point calculations. Additional methods available in Gaussian may be easily incorporated by simply updating gausslink with the relevant code.

Author Manuscript

The clear bottleneck for QM/MM/MC simulations is the QM calculations. When the previous BOSS–GAMESS interface was first developed, much effort was spent attempting to accelerate the code.[10] In that work and noted elsewhere,[14] significant speedup is achieved when a wavefunction from a previous QM calculation is used as an initial guess for the next QM calculation. This feature is implemented in the present interface by storing the Gaussian checkpoint file following each QM calculation and adding the Gaussian “GUESS=TCHECK” keyword to subsequent input files. The resultant reduction in computer time was approximately 10% for the present simulations of the Henry reaction in water. Throughput may be further increased by running the QM calculations in parallel on several CPUs.[14,17] When QM packages support parallel processing, the linker script may be edited to request use of these features. For the same simulation of a Henry reaction, using 2–4 processors for each QM calculation resulted in speedups of 30–35%. Together, both techniques delivered overall speedups of 40–45% for the Henry reaction. These measures are highly desirable for performing QM/MM simulations with the BOSS–Gaussian interface


Vilseck et al.

Page 6

Author Manuscript

in acceptable time lengths. Additional enhancements that have been pursued by others include use of approximate ab initio methods, such as resolution of the identity MP2 (RIMP2),[59] and messaging passing interface (MPI) assisted QM–MM data communication.[11,13,14,15] Monte Carlo/Free-Energy Perturbation Calculations (MC/FEP)

Author Manuscript

In BOSS, reaction modeling often involves merging the QM/MM/MC philosophy with Free Energy Perturbation (FEP) calculations.[4,5] FEP theory is employed to compute free energy changes along reaction coordinates, such as bonds that form or break during a reaction. The Zwanzig equation is used,[60] whereby free energy differences are calculated as the ensemble average of the potential energy difference between an initial (A) and final state (B) (eq. 3).[4,5,61–63] For reactions, inter- or intramolecular coordinates are perturbed.[4,47,51] To achieve sufficient convergence, the total perturbation is divided into a series of λ-windows, with geometrical changes of 0.01–0.05Å for bond length and 1–3° for bond angle perturbations per λ-window. Collectively, λ-windows are summed together to produce a free energy surface, or “potential of mean force” (PMF).[4] In the present work, both 1dimensional (1D) and 2-dimensional (2D) free energy surfaces have been generated from QM/MM/MC/FEP simulations, where either 1 or 2 degrees of freedom are independently perturbed.

(3)

Application to Three Reactions

Author Manuscript

Reaction modeling has a diverse history in developing QM/MM methods[1–4] and makes for an important testing ground to validate the BOSS–Gaussian interface. Experimental data or prior computational investigations are often available for comparisons. Three reactions were chosen to evaluate the performance of the BOSS–Gaussian interface. The Henry reaction between the anion of nitromethane and formaldehyde (HEN) was chosen for study first (Scheme 2A).[10,38] The small QM system size requirements, its 1D reaction profile, and prior investigation with AM1 and MP2 methods make it an appealing test. Two methyl transfer reactions between dimethylamine and dimethylammonium ion (MTA) and between dimethylamine and trimethylsulfonium ion (MTC) were next investigated (Scheme 2B and 2C).[64] In the prior PDDG/PM3 investigations of these reactions the computed activation barriers were over-estimated compared to experiment.[38,64]

Author Manuscript

Computational Details All Monte Carlo simulations, OPLS-AA energy evaluations, and FEP calculations were carried out in a modified version of BOSS 4.9.[65] Gaussian 09, revision.D.01, was used for the requisite QM single-point calculations for the solutes. All QM calculations employed the MP2 electronic structure method[55] and CM5 charges.[43] The QM calculations were run in parallel over 2–8 CPUs with at least 2 GB of memory.


Vilseck et al.

Page 7

Author Manuscript Author Manuscript

For all MC simulations, sampling was performed in the isothermal-isobaric ensemble at 25 °C and 1 atm. The aqueous environment was represented using 740 TIP4P water molecules with periodic boundary conditions.[48] Initial solvent coordinates originated from an equilibrated water box stored in BOSS with box dimensions of ca. 25 x 25 x 37.5 Å. All solute degrees of freedom were fully sampled except for the perturbed reaction coordinates. Solute moves were attempted every 100 configurations with molecular translations and rotations within ranges of ±0.06 Å and ±6.0°, respectively. TIP4P water molecules were freely translated and rotated within ranges of ±0.12 Å and ±12.0°, respectively, and all intramolecular degrees of freedom were fixed. The frequency of volumes moves was automatically adjusted by BOSS to yield average MC acceptance rates of 40%. Intermolecular cutoff distances were truncated at 12 Å, based on non-hydrogen atom pair distances, and interaction energies were quadratically smoothed to zero over the last 0.5 Å. All nonbonded interactions between QM and MM regions were calculated with the OPLSAA force field.[30,31] All solute-solvent Coulombic interactions employed unscaled CM1,[44] CM3,[45] or CM5[43] partial atomic charges, and standard OPLS-AA force field parameters were used for Lennard Jones interactions. Statistical uncertainties (± 1σ) were calculated using the batch means method,[61,66] with batch sizes of 500 K configurations. For the present reactions, standard uncertainties of 0.2–0.4 kcal/mol were obtained for free energies of activation and reaction. For each FEP window, 2.5–3.5 M configurations of equilibration and 5.0 M configurations of averaging were used (2.5–3.5 M/5.0 M). A cratic correction of 1.89 kcal/mol was added to each calculated absolute free energy (ΔG) to correspond to 1M standard states.[64,67,68]

The Henry Reaction Author Manuscript

As a fundamental C–C bond forming process, the Henry reaction is an interesting system to study.[69] Observed benefits for performing this reaction in water include short reaction times, facile purification, and reduced side-product formation.[70–72] Because the reaction is initiated by base-catalyzed proton abstraction from the nitroalkane,[73] most theoretical investigations have studied the reaction between the resultant nitroalkane anion and an aldehyde.[74,75] Although no experimental data exists for the parent reaction between nitromethane and formaldehyde, it is an appealing model due to the small size and straightforward reaction mechanism. The present investigation also mimics the validation performed for the BOSS-GAMESS interface.[10]

Author Manuscript

For the Henry reaction the forming C–C bond was chosen as the reaction coordinate.[10,38] Forty λ-windows with Δr(C–C) = 0.025 Å were used in conjunction with double-wide sampling (DWS)[61,62] to compute the total free energy profile for r(C–C) spanning 1.50 to 3.50 Å. This covers the reactant complex (RC), transition state (TS), and product (P) structures along the reaction path. Although AM1 SQM results were published in 2012,[38] AM1 calculations with CM1 charges (AM1/CM1A) were repeated here alongside MP2 calculations. The MP2/6-31+G(d) level of theory was used via the BOSS–Gaussian interface with CM5 charges (MP2/CM5).[55,76] This mirrors the level of theory used previously with the BOSS–GAMESS interface for the same reaction, only with CM2 charges (MP2/ CM2).[10,77] For efficiency, λ-windows were run in parallel on separate CPU nodes. Each window required ca. 150–175 K QM single-point calculations, totaling ca. 6–7 M QM


Vilseck et al.

Page 8

Author Manuscript

calculations for the full reaction profile. A single window of 2.5 M/5.0 M configurations took ca. 6.6 days to complete on an Intel Core 2 Quad 3.3 GHz processor.

Author Manuscript

Figure 1 shows the new free energy profiles for both AM1/CM1A and MP2/CM5 QM/MM/MC/FEP simulations. The free energies of activation (ΔG‡) are 15.3 and 11.3 kcal/mol for AM1/CM1A and MP2/CM5, respectively. With the BOSS–Gaussian interface, the MP2/CM5 transition state occurs at a C-C separation of 2.15 Å, which is earlier than the 1.95 Å with AM1/CM1A. This elongation was also observed with the BOSS–GAMESS calculations, but with a barrier height of 13.8 kcal/mol (see Table S1 and Figure S1).[10] These results can be compared with the optimized MP2/6-31+G(d) C–C separation of 2.29 Å for the gas-phase TS.[75] In both AM1/CM1A and MP2/CM5 simulations, the 2nitroethoxide product has r(C–C) = 1.55 Å; however, the predicted free energies of reaction (ΔGrxn) differ significantly. With AM1/CM1A, the reaction is predicted to be endergonic by 5.3 kcal/mol. With MP2/CM5, the reaction is predicted to be exergonic by −4.9 kcal/mol. With the BOSS–GAMESS calculations, a ΔGrxn of −4.7 kcal/mol was computed, matching well the present MP2 results.

Author Manuscript

The origin of the differences in free energies profiles warrants analysis. Reported gas-phase QM calculations using AM1, MP2, and CCSD(T) quantum methods offer insight into the AM1 versus MP2 difference.[38,74,75] Zorn et al. reported gas-phase CCSD(T)/aug-ccpVDZ//MP2/6-31+G(d) energies of activation and reaction of +3.2 and −1.4 kcal/mol. Zeropoint corrected energies were +2.3 and −2.6 kcal/mol, respectively.[75] These values match well the gas-phase free energies of activation and reaction by Kostal et al.[38] Relative to the reactant complex, they reported +3.9 and −0.5 kcal/mol, respectively, with MP2/6-31+G(d,p). The consistency of the MP2 and CCSD(T) results indicates that the reaction to form 2-nitroethoxide is exothermic in the gas-phase. Associated AM1 calculations by Kostal et al., however, produced gas-phase free energies of activation and reaction of +4.0 and +1.6 kcal/mol. The AM1 free energy of reaction is too positive; however, the 2.0 kcal/mol gas phase difference accounts for only a fraction of the total QM/MM ΔΔGrxn of 10.2 kcal/mol. The remaining difference must originate with the atomic charge distributions and differential solvation along the reaction path.

Author Manuscript

Differences in solvation can be assessed by perturbing from the CM1A to CM5 charges using a standard MC/FEP protocol.[33,46,47,51] This is a little complicated by the fact that the AM1/CM1A and MP2/CM5 simulations sample different solute geometries. By freezing all intramolecular motions of the solute and zeroing all intra-solute energetics, a calculation of the ΔΔGhyd(CM1A→CM5) should reflect only the CMx difference. In an attempt to represent the RC, TS, or P configurations that may be sampled in a full MC simulation, ten random starting configurations were used to compute ten ΔΔGhyd(CM1A→CM5) values for each RC, TS, and P stationary structure. Starting solute coordinates were generated every 2000 configurations from an AM1 QM/MM/MC gas-phase simulation. Solute geometries were constrained rigid and the CMx charges were obtained for each structure prior to the MC/FEP simulation. Charges were perturbed over 11 simple overlap sampling (SOS)[61,63] λ-windows in a box of 740 TIP4P water molecules. Computed ΔΔGhyd(CM1A→CM5) results were then averaged producing mean free energy differences of +5.69, +1.86, and −0.71 kcal/mol for RC, TS, and P structures, respectively. Thus, the RC and TS structures J Comput Chem. Author manuscript; available in PMC 2016 October 15.

Vilseck et al.

Page 9


are better solvated with CM1A charges and the product is a little better solvated with CM5 charges. The excessive endergonicity for the AM1/CM1A results stems primarily from too favorable hydration of the reactants and transition state, and not from poor solvation of the product. Figure 2 shows averaged CM1A and CM5 charges for all heavy atoms from the rigid-solute MC/FEP simulations. As expected, charge transfer is observed to flow from the nitro group of nitromethane anion to the carbonyl oxygen of formaldehyde as the reaction progresses from RC to P. On average, the CM1A charges for the reacting anion are more polarized by −0.3, +0.7, and −0.2 electrons than with CM5 charges for C, N, and O, respectively. This strengthens the interactions of the AM1/CM1A RC with surrounding water molecules. From the fixed solute MC/FEP simulations, a ΔΔG‡ and ΔΔGrxn(CM1A→CM5) of −3.8 and −6.4 kcal/mol can be approximated. This largely accounts for the ΔΔG‡ and ΔΔGrxn(AM1/CM1A→MP2/CM5) differences of −4.0 and −10.2 kcal/mol in Figure 1. Thus, different representations of partial atomic charges lead to the majority of the differences in the QM/MM reaction profiles in Figure 1.

Author Manuscript

Solute–solvent energy pair distribution functions (EPDs) can also be used to characterize the changes in solvation along the reaction path. EPDs record the average number of solute– solvent pair interactions and their respective energies. AM1/CM1A and MP2/CM5 EPDs are graphically displayed in Figure 3. The large band near zero represents the many weak interactions between solutes and distant bulk water molecules. Hydrogen bonding interactions appears to the left of this band with interaction energies between −20.0 and −5.0 kcal/mol. In accord with expected weakening of hydrogen bonds to the charge-delocalized TS, the RC band, located near −12.5 kcal/mol, shifts right by 3.0 kcal/mol to −9.5 kcal/mol as the reaction approaches the TS. Upon crossing the TS, the hydrogen-bonding interactions strengthen from −9.5 kcal/mol to ca. −15.0 kcal/mol for the charge localized 2-nitroethoxide product. Much of the negative charge is now concentrated on the ethoxide oxygen, producing highly favorable hydrogen bonds to water. Consistent with the differences in Figure 1, the hydrogen-bonding band for RC is broader with the AM1/CM1A calculations, while the product band is shifted to lower energy with MP2/CM5. In summary, the observed ΔΔGrxn of 10.2 kcal/mol is a consequence of both inherent AM1 versus MP2 method differences and the CM1A versus CM5 charges. However, the charge differences contribute to ΔΔGrxn to a greater extent for this reaction. Though no experimental reference exists to clarify whether the aqueous-phase reaction is endergonic or exergonic, PCM/MP2 studies by Kostal,[10] and addition of a single nitromethane solvent molecule into MP2 gas-phase calculations by Zorn et al.[75] both suggest that the reaction is exothermic in solution.

Author Manuscript

Methyl Transfer Reactions The transfer of a methyl group from ammonium or sulfonium ions to neutral amines follows an SN2 mechanism.[64] Further testing of the BOSS–Gaussian approach was performed for the two methyl transfer reactions. They are attractive and challenging because sufficient experimental and theoretical data are available for quantitative comparisons,[64,78] and 2DPMFs are required for the identification of the transition structures and activation barriers. Furthermore, previous PDDG/PM3 SQM investigations over-estimated the reaction barriers


Vilseck et al.

Page 10

Author Manuscript

in water by ca. 13.0 kcal/mol for three methyl transfer reactions, including MTA and MTC.[64]

Author Manuscript

For both reactions, the bond-forming and bond-breaking distances were defined as the two reaction coordinates, r1 and r2.[64] Specifically, N1–C and N2–C bonds, for MTA, and S–C and N–C bonds, for MTC, were chosen (see Scheme 2B and 2C). For each reaction, two 2DPMF surfaces were generated to locate the aqueous-phase transition states. The first surface is used to determine the general location of the TS, while the second pinpoints its exact position.[4] For both MTA and MTC, increments of Δr1 = 0.03 Å and Δr2 = 0.02 Å were used to generate the first 2D-PMFs. Bond lengths sampled were r1(N1–C) = 1.87 – 2.11 Å and r2(N2–C) = 1.96 – 2.12 Å for MTA, and r1(S–C) = 2.15 – 2.39 Å and r2(N–C) = 2.02 – 2.22 Å for MTC. The second 2D-PMFs utilized step sizes of Δr1 = Δr2 = 0.01 Å with bond length ranges of r1(N1–C) = r2(N2–C) = 1.90 – 2.00 Å for MTA, and r1(S–C) =2.16 – 2.26 Å and r2(N–C) = 2.08 – 2.12 Å for MTC. From the latter 2D-PMFs, TS geometries were identified as the highest point along the minimum free energy path (see Figure 4 and Figures S2–S4). Once a transition states was pinpointed, a 1D-PMF was calculated perturbing each TS to the corresponding RC, thus determining the free energy of activation.

Author Manuscript

Computational requirements are much larger for these reactions compared to the Henry reaction, owing to their increased size and the need for 2D-PMFs. MP2/6-31G(d) calculations were chosen as a reasonable balance between accuracy and speed.[55,76] Five to ten DWS windows were used for each 1D component of the full 2D surface. λ-Windows were run in parallel and required approximately 180 K QM single-point calculations per window, totaling ca. 20 M QM calculations for each reaction. On an Intel Core 2 Quad 3.3 GHz processor, a single window took between 6–10 days for 3.5 M/5.0 M configurations of sampling, compared to 10–15 minutes with PDDG/PM3. For completeness and consistency, PDDG/PM3 calculations with CM3P charges (PDDG/CM3) were recomputed along with the MP2/CM5 calculations. The repeated PDDG/CM3P results match well the results from 2007.[64] Differences of 0.05 Å or less in locations of transition states were found. For example, the MTA TS was previously located at r1 = 1.99 Å and r2 = 1.95 Å,[64] while the present calculations gave r1 = r2 = 2.00 Å. The MTC TS was identified at r1 = 2.28 Å and r2 = 2.04 Å. The newly computed ΔG‡s for MTA and MTC are 51.0 and 41.6 kcal/mol, respectively, which compare reasonably well with the previously published values of 49.2 and 39.3 kcal/mol.[64]

Author Manuscript

The MP2/CM5 QM/MM TS geometries are very similar to the PDDG/CM3P structures. As shown in Figure 4, the MTA TS was found at r1 = 1.96 Å and r2 = 1.97 Å. In view of the PDDG/CM3P and MP2/CM5 results, the MTA transition state appears to have nearly equal N–C bond lengths of 1.95 – 2.00 Å. For MTC, the total S–N separation is almost identical with both methods, ca. 4.3 Å, but the methyl group is farther from sulfur with PDDG/ CM3P. With MP2/CM5, the S–C bond length decreases from the PDDG/CM3P length of 2.28 Å to 2.19 Å. This is offset by an elongated N–C bond length of 2.10 Å, compared to the PDDG/CM3P value of 2.04 Å. Subsequent calculations perturbing the transition states to reactant complexes with r2 = 5.50 Å produced final free energies of activation of 29.2 and 21.3 kcal/mol for the MTA and MTC reactions. The experimental free energies of activation


Vilseck et al.

Page 11

Author Manuscript

are 34.4 and 28.1 kcal/mol, respectively.[78g] Agreement with experiment is improved using the MP2/CM5 calculations by ca. 6.0 kcal/mol, however, the computed barriers now underestimate the experimental heights.

Author Manuscript

As with the Henry reaction, the large ΔΔG‡ of ca. 19.0 kcal/mol must originate from the use of different CMx charge models and QM methods in the QM/MM simulations. Charge model differences are considered first. Tables S2 and S3 report final atomic charges for all MTA and MTC atoms taken from the last batch of MC sampling in the RC and TS λwindows; Tables S4 and S5 abbreviate this list by summing hydrogen charges into neighboring heavy atoms. In general, minor differences exist between the CM3P and CM5 charges. Given similar patterns of charge distribution, is it unsurprising that calculated energy pair distributions are similar between PDDG/CM3P and MP2/CM5 simulations (Figures S5 and S6). In progressing from reactants to the transition state, both PDDG/CM3P and MP2/CM5 simulations indicate a loss of strong solute–water interactions. Additional rigid solute MC/FEP calculations perturbing from CM3P to CM5 charges show only small differences between the CMx models. Following the same methodology as above for the Henry reaction, averaged free energies differences of 9.78 and 5.42 kcal/mol for MTA and −0.53 and 0.25 for MTC were computed for rigid RC and TS perturbations, respectively. These yield ΔΔG‡(CM3P→CM5) estimates of −4.36 and 0.78 kcal/mol for the MTA and MTC reactions. Collectively, these results indicate that the majority of the ca. 19 kcal/mol ΔΔG‡(PDDG/CM3P→MP2/CM5) is not a result of using different CMx charge models, but rather stems from the differences in the QM methodologies.

Author Manuscript

As part of the original 2007 publication,[64] CBS-QB3 computed gas-phase activation enthalpies for MTA and MTC were 9.9 and 7.9 kcal/mol; MP2/6-31G(d) enthalpies were 8.7 and 5.9 kcal/mol; and PDDG/PM3 enthalpies were 31.3 and 26.8 kcal/mol, respectively. On average, the PDDG/PM3 gas phase activation enthalpies were overestimated by ca. 20.0 kcal/mol, whereas MP2/6-31G(d) underestimated the barriers by 1.2–2.0 kcal/mol. These patterns of over- or under-estimation are maintained in the present QM/MM investigations. The PDDG/PM3 error is related to its too positive heats of formation for ammonium and sulfonium cations, which increases with branching (see Tables S6 and S7). For species involved in the MTA or MTC reactions, the error in heat of formation is largest for trimethylammonium ion, 18.1 kcal/mol. The error likely carries over to the two transition states leading to its formation.

Author Manuscript

The closer agreement between the MP2/CM5 results and experiment is encouraging. Table 1 compares the present MP2/CM5 QM/MM results with ab initio results using the CPCM and GB/SA continuum models from the 2007 paper.[64] The computed values are reasonably consistent. The principal anomaly is the erroneous prediction of the higher barrier for the MTC reaction than MTA with the GB/SA model. Improving the MP2/CM5 agreement with experiment may be possible by using larger basis sets. Gunaydin et al. report that diffuse functions are valuable for modelling methyl transfer reactions. For gas-phase reactions, they observed differences in activation enthalpies of 0.9– 2.6 kcal/mol between 6–31G(d) and 6–31+G(d,p) basis sets with B3LYP and MP2 methods.[64] Investigations of larger basis sets were not pursued here due to the increased


Vilseck et al.

Page 12

Author Manuscript

computational demands. For example, a single FEP window covering 2.5 M/5.0 M configurations of sampling with MP2/6-31+G(d,p)/CM5 is estimated to take more than 20 days on an Intel Xeon 2.66 GHz processor. In cases like these where semiempirical and ab initio transition structures for a reaction in the gas-phase are similar,[64] a cost-saving approach might begin with locating the TS in solution using SQM-based FEP calculations. Then, the TS location could be refined and final free energies of activation and reaction determined using more rigorous ab initio or density functional methods in the QM/MM simulations.

Conclusions

Author Manuscript Author Manuscript Author Manuscript

An interface between the BOSS and Gaussian programs has been developed. The BOSS program was modified to request quantum mechanical energies and partial atomic charges from Gaussian 09. In this manner, CM5 charges and ab initio or density functional methods are incorporated into BOSS QM/MM capabilities. Data communication between the programs is controlled by an external C-shell linker script. Given sufficient computational resources, any option for a single-point calculation in Gaussian is available for use in the interface. However, as millions of QM calculations are normally required for a full QM/MM/MC/FEP study, applied levels of theory need to consider a balance between accuracy and speed. Simulation accelerations of ca. 40% were observed by reading in a previous wave function as an initial guess and running QM calculations in parallel. Nevertheless, in the present study, a single FEP λ-window took on average between 6–10 days with 3.5 M/5.0 M configurations of MC sampling and parallel processing. In comparison, AM1 and PDDG/PM3 semiempirical QM methods required less than 20 minutes per FEP λ-window on a single processor. To validate the operation and application of the interface, three organic reactions in water were examined. Significant differences between MP2 and SQM computed free energies of activation and reaction were observed. Differences in CMx charges largely affected computed free energies of reaction for the Henry reaction. For the methyl transfer reactions, however, free energy discrepancies were mainly a result of differences in the intrinsic MP2 and PDDG/PM3 results. Inclusion of ab initio or density function methods and CM5 charges in on-the-fly QM/MM simulations provides the opportunity for generating much more accurate results for a wide range of problems including reaction profiles. Cm5 charges performed well and are recommended for continued QM/MM use. Finally, the similarity between MP2 and PDDG/PM3 transition states for the MTA and MTC reactions suggests that the more cost-effective application of the BOSS–Gaussian interface is in computing the final free-energy profile, not in the initial two-dimensional scan for locating the TS. When gas-phase ab initio and SQM TS geometries are similar, SQM methods could be used for initial mapping of the free-energy surface in solution. The ab initio or DFT-based QM/MM calculations can then be used to refine the TS location and to determine the free energies of activation and reaction.

Supplementary Material Refer to Web version on PubMed Central for supplementary material.


Vilseck et al.

Page 13

Author Manuscript

Acknowledgments Gratitude is expressed to the National Institutes of Health (GM32136) for support. This work was also made possible by the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center. Gratitude is expressed to John Faver, Daniel Cole, Michael Robinson, and Leela Dodda for helpful discussions. This paper is dedicated to the memory of Professor Paul von Ragué Schleyer, an inspirational teacher and extraordinarily creative scientist.

References


1. (a) Warshel A, Karplus M. J Am Chem Soc. 1972; 94:5612–5625.(b) Warshel A, Levitt M. J Mol Biol. 1976; 103:227–249. [PubMed: 985660] 2. (a) Karplus M. Angew Chem Int Ed. 2014; 53:9992–10005.(b) Levitt M. Angew Chem Int Ed. 2014; 53:10006–10018.(c) Warshel A. Angew Chem Int Ed. 2014; 53:100120–10031.(d) Jorgensen WL. Cell. 2013; 155:1199–1202. [PubMed: 24315087] 3. Some review articles: Acevedo O. J Phys Chem A. 2014; 118:11653–11666. [PubMed: 25329366] Liu M, Wang Y, Chen Y, Field MJ, Gao J. Isr J Chem. 2014; 54:1250–1263.van der Kamp MW, Mulholland AJ. Biochemistry. 2013; 52:2708–2728. [PubMed: 23557014] Steinbrecher T, Elstner M. Biomolecular Simulations. Monticelli L, Salonen E. Methods in Molecular Biology. Humana PressNew York2013:91–124.Lodola A, De Vivo M. Adv Protein Chem Struct Biol. 2012; 87:337– 362. [PubMed: 22607760] Sousa SF, Fernandes PA, Ramos MJ. Phys Chem Chem Phys. 2012; 14:12431–12441. [PubMed: 22870506] Senn HM, Thiel W. Angew Chem Int Ed. 2009; 48:1198– 1229.Hu H, Yang W. Annu Rev Phys Chem. 2008; 59:573–601. [PubMed: 18393679] Lin H, Truhlar DG. Theor Chem Acc. 2007; 117:185–199.Friesner RA, Guallar V. Annu Rev Phys Chem. 2005; 56:389–427. [PubMed: 15796706] Gao J, Truhlar DG. Annu Rev Phys Chem. 2002; 53:467– 505. [PubMed: 11972016] Gordon MS, Schmidt MW. Computational Science – ICCS 2003. Sloot PMA, Abramson D, Bogdanov AV, Gorbachev YE, Dongarra JJ, Zomaya AY. Lecture Notes in Computer Science. SpringerBerlin Heidelberg2003:75–83. 4. (a) Acevedo O, Jorgensen WL. WIREs Comput Mol Sci. 2014; 4:422–435.(b) Acevedo O, Jorgensen WL. Acc Chem Res. 2010; 43:142–151. [PubMed: 19728702] 5. Jorgensen WL, Tirado-Rives J. J Comput Chem. 2005; 26:1689–1700. [PubMed: 16200637] 6. Dewar MJS, Zoebisch EG, Healy EF, Stewart JJP. J Am Chem Soc. 1985; 107:3902–3909. 7. (a) Stewart JJP. J Comput Chem. 1989; 10:209–220.(b) Stewart JJP. J Comput Chem. 1989; 10:221–264. 8. (a) Repasky MP, Chandrasekhar J, Jorgensen WL. J Comput Chem. 2002; 23:1601–1622. [PubMed: 12395428] (b) Tubert-Brohman I, Guimaraes CRW, Repasky MP, Jorgensen WL. J Comput Chem. 2004; 25:138–150. [PubMed: 14635001] (c) Tubert-Brohman I, Guimaraes CRW, Jorgensen WL. J Chem Theory Comput. 2005; 1:817–823. [PubMed: 19011692] 9. (a) Field MJ, Bash PA, Karplus M. J Comput Chem. 1990; 11:700–733.(b) Lyne P, Hodoscek M, Karplus M. J Phys Chem A. 1999; 103:3462–3471.(c) Thiel W, Voityuk AA. J Phys Chem. 1996; 100:616–626.(d) Luque FJ, Reuter N, Cartier A, Ruiz-Lopez MF. J Phys Chem A. 2000; 104:10923–10931.(e) Stewart JJP. J Computer-Aided Mol Design. 1990; 4:1–105.(f) Rzepa HS, Yi M. J Chem Soc, Perkin Trans 2. 1990:943–951.(g) Jurema MW, Shields GC. J Comput Chem. 1993; 14:89–104. 10. Kostal, J. Ph D Thesis. Yale University; New Haven, CT: 2012. 11. Riahi S, Rowley CN. J Comput Chem. 2014; 35:2076–2089. [PubMed: 25178266] 12. Woodcock HL, Hodoscek M, Gilbert ATB, Gill PMW, Schaefer HF, Brooks BR. J Comput Chem. 2007; 28:1485–1502. [PubMed: 17334987] 13. Loferer MJ, Loeffler HH, Liedl KR. J Comput Chem. 2003; 24:1240–1249. [PubMed: 12820132] 14. (a) Gotz AW, Clark MA, Walker RC. J Comput Chem. 2014; 35:95–108. [PubMed: 24122798] (b) Isborn CM, Gotz AW, Clark MA, Walker RC, Martinez TJ. J Chem Theory Comput. 2012; 8:5092–5106. [PubMed: 23476156] 15. Meier K, Schmid N, van Gunsteren WF. J Comput Chem. 2012; 33:2108–2117. [PubMed: 22736402]


Vilseck et al.

Page 14

Author Manuscript Author Manuscript Author Manuscript Author Manuscript

16. (a) Torras J, Seabra GM, Deumens E, Trickey SB, Roitberg AE. J Comput Chem. 2008; 29:1564– 1573. [PubMed: 18270957] (b) Roberts BP, Seabra SM, Roitberg AE, Merz KM, Deumens E, Torras J, Trickey SB. J Comput Chem. 2013; 33:1643–1644. [PubMed: 22570199] (c) Torras, J.; Deumens, E.; Trickey, SB.; Cheng, HP.; Cao, C.; He, Y.; Muralidharan, K.; Roitberg, A.; Seabra, G.; Roberts, BP.; Bertran, O.; Fu, Z. PUPIL, Program for User Package Interface and Linking, version 3.0. Universitat Politècnica de Catalunya; Spain: 2014. 17. Okamoto T, Yamada K, Koyano Y, Asada T, Koga N, Nagaoka M. J Comput Chem. 2011; 32:932–942. [PubMed: 20949515] 18. Hagiwara Y, Ohta T, Tateno M. J Phys: Condens Matter. 2009; 21:064234. [PubMed: 21715936] 19. Biswas PK, Gogonea V. J Chem Phys. 2005; 123:164114. [PubMed: 16268688] 20. (a) Metz S, Kastner J, Sokol AA, Keal TW, Sherwood P. WIREs Comput Mol Sci. 2014; 4:101– 110.(b) Sherwood P, de Vries AH, Guest MF, Schreckenbach G, Catlow CRA, French SA, Sokol AA, Bromley ST, Thiel W, Turner AJ, Billeter S, Terstegen F, Thiel S, Kendrick J, Rogers SC, Casci J, Watson M, King F, Karlsen E, Sjovoll M, Fahmi A, Schafer A, Lennartz C. J Mol Struct: THEOCHEM. 2003; 632:1–28. 21. (a) Wang J, Wolf RM, Caldwell JW, Kollman PA, Case DA. J Comput Chem. 2004; 25:1157– 1174. [PubMed: 15116359] (b) Salomon-Ferrer R, Case DA, Walker RC. WIREs Comput Mol Sci. 2013; 3:198–210. 22. (a) Scott WRP, Hunenberger PH, Tironi IG, Mark AE, Billeter SR, Fennen J, Torda AE, Huber T, Kruger P, van Gunsteren WF. J Phys Chem A. 1999; 103:3596–3607.(b) Christen M, Hunenberger PH, Bakowies D, Baron R, Burgi R, Geerke DP, Heinz TN, Kastenholz MA, Krautler V, Oostenbrink C, Peter C, Trzesniak D, van Gunsteren WF. J Comput Chem. 2005; 26:1719–1751. [PubMed: 16211540] 23. (a) Brooks BR, Brooks CL, Mackerell AD, Nilsson L, Petrella RJ, Roux B, Won Y, Archontis G, Bartels C, Boresch S, Caflisch A, Caves L, Cui Q, Dinner AR, Feig M, Fischer S, Gao J, Hodoscek M, Im W, Kuczera K, Lazaridis T, Ma J, Ovchinnikov V, Paci E, Pastor RW, Post CB, Pu JZ, Schaefer M, Tidor B, Venable RM, Woodcock HL, Wu X, Yang W, York DM, Karplus M. J Comput Chem. 2009; 30:1545–1614. [PubMed: 19444816] (b) Vanommeslaeghe K, Hatcher E, Acharya C, Kundu S, Zhong S, Shim J, Darian E, Guvench O, Lopes P, Vorobyov I, Mackerell AD. J Comput Chem. 2010; 31:671–690. [PubMed: 19575467] (c) Momany FA, Rone R. J Comput Chem. 1992; 13:888–900.(d) Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. J Comput Chem. 1983; 4:187–217. 24. Ponder, JW. TINKER – Software Tools for Molecular Design. Washington University School of Medicine; Saint Louis, MO: 2014. 25. Schmidt MW, Baldridge KK, Boatz JA, Elbert ST, Gordon MS, Jensen JH, Koseki S, Matsunaga N, Nguyen KA, Su SJ, Windus TL, Dupuis M, Montgomery JA. J Comput Chem. 1993; 14:1347– 1367. 26. Frisch, MJ.; Trucks, GW.; Schlegel, HB.; Scuseria, GE.; Robb, MA.; Cheeseman, JR.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, GA.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, HP.; Izmaylov, AF.; Bloino, J.; Zheng, G.; Sonnenberg, JL.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, JA.; Peralta, JE.; Ogliaro, F.; Bearpark, M.; Heyd, JJ.; Brothers, E.; Kudin, KN.; Staroverov, VN.; Keith, T.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, JC.; Iyengar, SS.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, NJ.; Klene, M.; Knox, JE.; Cross, JB.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, RE.; Yazyev, O.; Austin, AJ.; Cammi, R.; Pomelli, C.; Ochterski, JW.; Martin, RL.; Morokuma, K.; Zakrzewski, VG.; Voth, GA.; Salvador, P.; Dannenberg, JJ.; Dapprich, S.; Daniels, AD.; Farkas, Ö.; Foresman, JB.; Ortiz, JV.; Cioslowski, J.; Fox, DJ. Gaussian 09, revision D.01. Gaussian, Inc; Wallingford, CT: 2013. 27. (a) Ahlrichs R, Bar M, Haser M, Horn H, Kolmel C. Chem Phys Lett. 1989; 162:165–169.(b) Treutler O, Ahlrichs R. J Chem Phys. 1995; 102:346–354.(c) Furche F, Ahlrichs R, Hattig C, Klopper W, Sierka M, Weigend F. WIREs Comput Mol Sci. 2013; 4:91–100. 28. (a) Shao Y, Fusti-Molnar L, Jung Y, Kussmann J, Ochsenfeld C, Brown ST, Gilbert ATB, Slipchenko LV, Levchenko SV, O’Neill DP, Distasio RA, Lochan RC, Wang T, Beran GJO, Besley NA, Herbert JM, Lin CY, Van Voorhis T, Chien SH, Sodt A, Steele RP, Rassolov VA, Maslen PE, Korambath PP, Adamson RD, Austin B, Baker J, Byrd EFC, Daschel H, Doerksen RJ,


Vilseck et al.

Page 15


Dreuw A, Dunietz BD, Dutoi AD, Furlani TR, Gwaltney SR, Heyden A, Hirata S, Hsu CP, Kedziora G, Khalliulin RZ, Klunzinger P, Lee AM, Lee MS, Liang W, Lotan I, Nair N, Peters B, Proynov EI, Pieniazek PA, Rhee YM, Ritchie J, Rosta E, Sherrill DC, Simmonett AC, Subotnik JE, Woodcock HL, Zhang W, Bell AT, Chakraborty AK, Chipman DM, Keil FJ, Warshel A, Hehre WJ, Schaefer HF, Kong J, Krylov AI, Gill PMW, Head-Gordon M. Phys Chem Chem Phys. 2006; 8:3172. [PubMed: 16902710] (b) Krylov AI, Gill PMW. WIREs Comput Mol Sci. 2013; 3:317–326. 29. (a) Ufimtsev IS, Martinez TJ. J Chem Theory Comput. 2009; 5:1004–1015.(b) Ufimtsev IS, Martinez TJ. J Chem Theory Comput. 2009; 5:2619–2628. 30. Jorgensen WL, Tirado-Rives J. Proc Nat Acad Sci USA. 2005; 102:6665–6670. [PubMed: 15870211] 31. Jorgensen WL, Maxwell DS, Tirado-Rives J. J Am Chem Soc. 1996; 118:11225–11236. 32. Vilseck JZ, Sambasivarao SV, Acevedo O. J Comput Chem. 2011; 32:2836–2842. [PubMed: 21732390] 33. (a) Vilseck JZ, Tirado-Rives J, Jorgensen WL. J Chem Theory Comput. 2014; 10:2802–2812. [PubMed: 25061445] (b) Dodda LS, Vilseck JZ, Cutrona KJ, Jorgensen WL. 2015 submitted. 34. Vilseck JZ, Tirado-Rives J, Jorgensen WL. Phys Chem Chem Phys. 2015; 17:8407–8415. [PubMed: 25589343] 35. Wolters LP, Schyman P, Pavan MJ, Jorgensen WJ, Bickelhaupt FM, Kozuch S. WIREs Comput Mol Sci. 2014; 4:523–540. 36. Yan XC, Schyman P, Jorgensen WL. J Phys Chem A. 2014; 118:2820–2826. [PubMed: 24678636] 37. Jorgensen WL, Schyman P. J Chem Theory Comput. 2012; 8:3895–3901. [PubMed: 23329896] 38. Kostal J, Voutchkova AM, Jorgensen WL. Org Lett. 2012; 14:260–263. [PubMed: 22168236] 39. Alexandrova A, Jorgensen WL. J Phys Chem B. 2011; 115:13624–13632. [PubMed: 21995727] 40. Kostal J, Jorgensen WL. J Am Chem Soc. 2010; 132:8766–8773. [PubMed: 20524660] 41. Acevedo O, Jorgensen WL. J Phys Chem B. 2010; 114:8425–8430. [PubMed: 20527873] 42. Thomas LL, Tirado-Rives J, Jorgensen WL. J Am Chem Soc. 2010; 132:3097–3104. [PubMed: 20148559] 43. Marenich AV, Jerome SV, Cramer CJ, Truhlar DG. J Chem Theory Comput. 2012; 8:527–541. 44. (a) Storer JW, Giesen DJ, Cramer CJ, Truhlar DG. J Comput-Aided Mol Des. 1995; 9:87–110. [PubMed: 7751872] (b) Mulliken RS. J Chem Phys. 1935; 3:564–572.(c) Mulliken RS. J Chem Phys. 1955; 23:1833–1840.(d) Mulliken RS. J Chem Phys. 1962; 36:3428–3439. 45. (a) Thompson JD, Cramer CJ, Truhlar DG. J Comput Chem. 2003; 24:1291–1304. [PubMed: 12827670] (b) Löwdin PO. J Chem Phys. 1950; 18:365–375.(c) Baker J. Theor Chim Acta. 1985; 68:221–229. 46. Udier-Blagovic M, Morales De Tirado P, Pearlman SA, Jorgensen WL. J Comput Chem. 2004; 108:16264–16270. 47. Kaminski G, Jorgensen WL. J Phys Chem B. 1998; 102:1787–1796. 48. Jorgensen WL, Chandrasekhar J, Madura JD, Impey RW, Klein ML. J Chem Phys. 1983; 79:926– 935. 49. Metropolis N, Rosenbluth AW, Rosenbluth MN, Teller AH, Teller E. J Chem Phys. 1953; 21:1087–1092. 50. Guimaraes CRW, Udier-Blagovic M, Jorgensen WL. J Am Chem Soc. 2005; 127:3577–3588. [PubMed: 15755179] 51. Carlson HA, Nguyen TB, Orozco M, Jorgensen WL. J Comput Chem. 1993; 14:1240–1249. 52. (a) Roothaan CC. J Rev Mod Phys. 1951; 23:69–89.(b) Pople JA, Nesbet RK. J Chem Phys. 1954; 22:571–572. 53. (a) Becke AD. Phys Rev A. 1988; 38:3098–3100. [PubMed: 9900728] (b) Lee C, Yang W, Parr RG. Phys Rev B. 1988; 37:785–789.(c) Becke AD. J Chem Phys. 1993; 98:5648–5652.(d) Stephens PJ, Devlin FJ, Chabalowski CF, Frisch MJ. J Phys Chem. 1994; 98:11623–11627. 54. (a) Zhao Y, Truhlar DG. J Chem Phys. 2006; 125:194101. [PubMed: 17129083] (b) Zhao Y, Truhlar DG. Acc Chem Res. 2008; 41:157–167. [PubMed: 18186612] (c) Zhao Y, Truhlar DG. Theor Chem Acc. 2008; 120:215–241. J Comput Chem. Author manuscript; available in PMC 2016 October 15.

Vilseck et al.

Page 16


55. (a) Moller C, Plesset MS. Phys Rev. 1934; 46:618–622.(b) Head-Gordon M, Pople JA, Frisch MJ. Chem Phys Lett. 1988; 153:503–506.(c) Saebo S, Almlof J. Chem Phys Lett. 1989; 154:83–89.(d) Frisch MJ, Head-Gordon M, Pople JA. Chem Phys Lett. 1990; 166:275–280.(e) Frisch MJ, HeadGordon M, Pople JA. Chem Phys Lett. 1990; 166:281–289.(f) Head-Gordon M, Head-Gordon T. Chem Phys Lett. 1994; 220:122–128. 56. (a) Hirshfeld FL. Theor Chim Acta. 1977; 44:129–138.(b) Ritchie JP. J Am Chem Soc. 1985; 107:1829–1837.(c) Ritchie JP, Bachrach SM. J Comput Chem. 1987; 8:499–509. 57. (a) Singh UC, Kollman PA. J Comp Chem. 1984; 5:129–145.(b) Besler BH, Merz KM, Kollman PA. J Comp Chem. 1990; 11:431–439. 58. Breneman CM, Wiberg KB. J Comp Chem. 1990; 11:361–373. 59. (a) Feyereisen M, Fitzgerald G, Komornicki A. Chem Phys Lett. 1993; 208:359–363.(b) Vahtras O, Almlof J, Feyereisen MW. Chem Phys Lett. 1993; 213:514–518.(c) Bernholdt DE, Harrison RJ. Chem Phys Lett. 1996; 250:470–484. 60. Zwanzig RW. J Chem Phys. 1954; 22:1420–1426. 61. Jorgensen WL, Thomas LL. J Chem Theory Comput. 2008; 4:869–876. [PubMed: 19936324] 62. Jorgensen WL, Ravimohan C. J Chem Phys. 1985; 83:3050–3056. 63. Lu N, Kofke DA, Woolf TB. J Comput Chem. 2004; 25:28–39. [PubMed: 14634991] 64. Gunaydin H, Acevedo O, Jorgensen WL, Houk KN. J Chem Theory Comput. 2007; 3:1028–1035. 65. Jorgensen, WL. BOSS, version 4.9. Yale University; New Haven, CT: 2014. 66. Allen, MP.; Tildesley, DJ. Computer Simulations of Liquids. Clarendon; Oxford: 1987. p. 191-198. 67. (a) Yu BY. J Pharm Sci. 2001; 90:2099–2102. [PubMed: 11745768] (b) Hermans J, Wang L. J Am Chem Soc. 1997; 119:2707–2714.(c) Murphy KP, Xie D, Thompson KS, Amzel LM, Freire E. Proteins: Struct, Funct, Genet. 1994; 18:63–67. [PubMed: 8146122] 68. Vayner G, Houk KN, Jorgensen WL, Brauman JI. J Am Chem Soc. 2004; 126:9054–9058. [PubMed: 15264838] 69. Henry L, Seances CRH. Acad Sci. 1895; 120:1265–1267. 70. Ballini R, Bosica G. J Org Chem. 1997; 62:425–427. [PubMed: 11671420] 71. Phukan M, Borah KJ, Borah R. Synth Commun. 2008; 38:3068–3073. 72. Wang Z, Xue H, Wang S, Yuan C. Chem Biodiversity. 2005; 2:1195–1199. 73. Smith, MB.; March, J. March’s Advance Organic Chemistry: Reactions, Mechanisms, and Structure. 6. John Wiley & Sons, Inc; Hoboken, New Jersey: 2007. p. 1357-1358. 74. Lecea B, Arrieta A, Morao I, Cossio FP. Chem Eur J. 1997; 3:20–28. 75. Zorn D, Lin VSY, Pruski M, Gordon MS. J Phys Chem A. 2008; 112:10635–10649. [PubMed: 18823100] 76. (a) Ditchfield R, Hehre WJ, Pople JA. J Chem Phys. 1971; 54:724–728.(b) Hehre WJ, Ditchfield R, Pople JA. J Chem Phys. 1972; 56:2257–2261.(c) Hariharan PC, Pople JA. Theor Chem Acc. 1973; 28:213–22.(d) Hariharan PC, Pople JA. Mol Phys. 1974; 27:209–14.(e) Francl MM, Pietro WJ, Hehre WJ, Binkley JS, DeFrees DJ, Pople JA. J Chem Phys. 1982; 77:3654–3665.(f) Clark T, Chandrasekhar J, Spitznagel GW, Schleyer PvR. J Comput Chem. 1983; 4:294–301.(g) Frisch MJ, Pople JA, Binkley JS. J Chem Phys. 1984; 80:3265–3269. 77. Li J, Zhu T, Cramer CJ, Truhlar DG. J Phys Chem A. 1998; 102:1820–1831. 78. (a) Roca M, Marti S, Andres J, Moliner V, Tunon I, Bertran J, Williams IH. J Am Chem Soc. 2003; 125:7726–7737. [PubMed: 12812514] (b) Kuhn B, Kollman PA. J Am Chem Soc. 2000; 122:2586–2596.(c) Zheng YJ, Bruice TC. J Am Chem Soc. 1997; 119:8137–8145.(d) Lau EY, Bruice TC. J Am Chem Soc. 2000; 122:7165–7171.(e) Kahn K, Bruice TC. J Am Chem Soc. 2000; 122:46–51.(f) Castejon H, Wiberg KB, Sklenak S, Hinz W. J Am Chem Soc. 2001; 123:6092–6097. [PubMed: 11414843] (g) Callahan BP, Wolfenden R. J Am Chem Soc. 2003; 125:310–311. [PubMed: 12517124] 79. Leach S. Chem Phys. 2012; 392:170–179. 80. Chyall LJ, Byrd MHC, Kenttamaa HI. J Am Chem Soc. 1994; 116:10767–10772.66.


Vilseck et al.

Page 17


Figure 1.

AM1/CM1A (blue) and MP2/CM5 (red) QM/MM free energy reaction profiles for the Henry reaction between nitromethane anion and formaldehyde. The dashed line represents 0.00 kcal/mol; RC, TS, and P snapshots are illustrated.

Author Manuscript Author Manuscript J Comput Chem. Author manuscript; available in PMC 2016 October 15.

Vilseck et al.

Page 18


Figure 2.

Average AM1/CM1A and MP2/CM5 partial atomic charges for heavy atoms in the RC (normal), TS (bold), and P (italic) stationary structures for the Henry reaction obtained from ten rigid solute conformations.


Vilseck et al.

Page 19

Author Manuscript Figure 3.

Author Manuscript

Solute–water energy pair distribution functions for the Henry reaction between nitromethane anion and formaldehyde. AM1/CM1A results are on the left; MP2/CM5 is on the right. EPDs are colored by structure: RC (red), TS (blue), and P (green).


Vilseck et al.

Page 20

Author Manuscript Author Manuscript Figure 4.

Author Manuscript

The final MP2/CM5 2D-PMF for the MTA reaction. Relative free energies (kcal/mol) are colored from red (positive) to blue (negative). The star indicates the TS geometry.

Author Manuscript J Comput Chem. Author manuscript; available in PMC 2016 October 15.

Vilseck et al.

Page 21


Scheme 1.

Cyclic Workflow of the BOSS–Gaussian Interface. The programs that perform each operation are designated by color: BOSS (blue), Linker (green), and Gaussian (orange).

Author Manuscript J Comput Chem. Author manuscript; available in PMC 2016 October 15.

Vilseck et al.

Page 22

Author Manuscript Scheme 2.

Henry and Methyl Transfer Reactions Studied.

Author Manuscript Author Manuscript Author Manuscript J Comput Chem. Author manuscript; available in PMC 2016 October 15.

Vilseck et al.

Page 23

Table 1

Author Manuscript

Activation Barriers in Water for Two Methyl Transfer Reactions with QM/MM, CPCM, and GB/SA Solvation Models (kcal/mol)[a] Theory

Solvation Model

MTA

MTC

QM/MM/MC

29.2

21.3

CPCM

36.8

22.2

CBS-QB3

CPCM

36.3

23.8

CBS-QB3

GB/SA

29.9

34.6

34.4

28.1

MP2/6-31G(d)/CM5 MP2/6-31+G(d,p)

Experiment [a]

CPCM, GB/SA, and experimental data from Ref. 64.

Author Manuscript Author Manuscript Author Manuscript J Comput Chem. Author manuscript; available in PMC 2016 October 15.

MM and MD Simulations.

Computational analysis of methyl transfer reactions in dengue virus methyltransferase.

MM Molecular Dynamics Simulations of Enzyme Reactions.

Preparation of indium nitronates and their Henry reactions.

Computer simulations of electron-transfer reactions in solution and in photosynthetic reaction centers.

MM free energy simulations: recent progress and challenges.

MM MD simulations of La(III)-phosphopeptide complexes.

MM simulations of silica.

MM MD and free energy simulation study of methyl transfer processes catalyzed by PKMTs and PRMTs.

Trifluoromethyl-modified dipeptides by ZrCl4-promoted aza-Henry reactions.

MM simulations.

MM simulations.

MM Sampling of Enzymatic Reactions.

syn- and enantioselective Henry reactions of aliphatic aldehydes and application to the synthesis of safingol.

MM and MD Simulations.

MM methodology.

QM-MM simulations on p53-DNA complex: a study of hot spot and rescue mutants.

MM simulations of the obelin bioluminescence and fluorescence reveal an unexpected light emitter.

Gold-catalyzed carbene transfer reactions.

Layer Transfer from Chemically Etched 150 mm Porous Si Substrates.

Mesoscale simulations of shockwave energy dissipation via chemical reactions.

Investigating DNA Binding and Conformational Variation in Temperature Sensitive p53 Cancer Mutants Using QM-MM Simulations.

Pathway analysis of protein electron-transfer reactions.

(2S,5R)-2-Methylaminomethyl-1-methyl-5-phenylpyrrolidine, a chiral diamine ligand for copper(II)-catalysed Henry reactions with superb enantiocontrol.