Exo Selectivity of the Diels–Alder

May 29, 2019 - PDF (1 MB) .... reaction and three chemical reactions in aqueous solution. ...... Rasmussen, C.; Williams, C. Gaussian Processes for Ma...
4 downloads 0 Views 1MB Size
Download PDF
Article Cite This: J. Phys. Chem. B 2019, 123, 5131−5138

pubs.acs.org/JPCB

Computational Insights into Endo/Exo Selectivity of the Diels−Alder Reaction in Explicit Solvent at Ab Initio Quantum Mechanical/ Molecular Mechanical Level Pengfei Li,† Fengjiao Liu,‡ Yihan Shao,§ and Ye Mei*,†,∥,⊥,§

Downloaded via UNIV FRANKFURT on July 25, 2019 at 06:14:58 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

State Key Laboratory of Precision Spectroscopy, School of Physics and Materials Science, East China Normal University, Shanghai 200062, China ‡ College of Chemistry and Chemical Engineering, Shanghai University of Engineering Science, Shanghai 201620, China § Department of Chemistry and Biochemistry, University of Oklahoma, Norman, Oklahoma 73019, United States ∥ NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai 200062, China ⊥ Collaborative Innovation Center of Extreme Optics, Shanxi University, Taiyuan, Shanxi 030006, China ABSTRACT: The microscopic insight into the endo/exo stereoselectivity of the Diels−Alder (DA) reaction between cyclopentadiene and methyl vinyl ketone (MVK) has posed a challenge to the computational chemists, which requires an accurate free-energy (FE) landscape. Although ab initio (ai) quantum mechanical/molecular mechanical (QM/MM) calculations are capable of providing accurate conformation energies, they are far too expensive for free-energy calculations, in which large-scale sampling in the phase space is always indispensable. Recently, on the basis of the idea of reference-potential methods, a new approach termed MBAR+wTP method was proposed by us (Li; et al. J. Chem. Theory Comput. 2018, 14, 5583) for the calculations of ai QM/MM FE profiles with a much less computational expense. In this work, we applied this method to investigate the endo/exo stereoselectivity of the Diels−Alder (DA) reaction at B3LYP/MM level and the results indicate that this method can yield more accurate activation free energies than the semiempirical Hamiltonian PM6 and the “sampling-free” quantum mechanical methods at B3LYP and MP2 levels using a polarizable continuous solvent model (PCM). The stereoselectivity mainly comes from the solvation effect. At the transition state, the first peak of the solvent distribution near the oxygen atom in MVK is slightly closer for the endo pathway than that for the exo pathway, whereas the difference at PM6 level is negligible. Nonetheless, in terms of the stereoselectivity, we obtained an endo/exo ratio of 2 with this method, which is still 1 order of magnitude smaller than the experimental measurement. Although the endo/exo ratio from the MP2/PCM calculation agrees well with the experiment, the absolute reaction barriers deviate by about 5 kcal mol−1. This work issues a warning to the use of the computational methods that are routinely employed for the study of organic reactions in the condensed phase, and it also brings a new method that is affordable and more reliable.

■

10−5 M−1 s−1 and it decreases to 75.5 × 10−5 M−1 s−1 in methanol, corresponding to a difference in the reaction barrier of about 2 kcal mol−1.10,23 Moreover, this DA reaction prefers the endo isomer with a reaction rate of about 4400 × 10−5 M−1 s−1, compared with the exo isomer with a reaction rate of about 180 × 10−5 M−1 s−1,10,12,13,15 which results in a reaction rate ratio of about 25 for the endo and exo isomers. This small difference in the reaction barrier poses an enormous challenge to computational investigations. Jorgensen and co-workers found that the endo-cis (with MVK in an s-cis conformation) transition state is favored in the gas phase, where a reaction barrier of 29.42 kcal mol−1 for the endo-cis conformation is lower than that of 29.86 kcal mol−1 for the

INTRODUCTION The Diels−Alder (DA) reaction is particularly useful for forming unsaturated 6-membered ring systems in synthetic organic chemistry.1−4 It has been found that this class of reactions has not only, in general, concerted mechanisms5,6 but may also follow stepwise mechanisms.7−9 The Diels−Alder reaction between cyclopentadiene (CP) and methyl vinyl ketone (MVK), in which CP is a conjugated diene and MVK is a dienophile, as shown in Figure 1, was first reported by Breslow et al. in 1980.10 Wide studies from experimental10−16 and computational aspects17−22 suggest that the acceleration of DA reactions in solutions stems from enhanced polarization of the transition state (TS), which leads to stronger hydrogen bonding between the solvent and the TS structure. Breslow and co-workers found that not only the reactivity but also the endo/exo selectivity in this DA reaction can be enhanced by water.10,13,14 The reaction rate in water is 4400 × © 2019 American Chemical Society

Received: March 1, 2019 Revised: April 29, 2019 Published: May 29, 2019 5131

DOI: 10.1021/acs.jpcb.9b01989 J. Phys. Chem. B 2019, 123, 5131−5138

Article

The Journal of Physical Chemistry B

Figure 1. Diels−Alder reaction between cyclopentadiene (CP) and methyl vinyl ketone (MVK), where η1 = dC1C3 and η2 = dC2C4 were chosen as the two-dimensional (2D) reaction coordinates and both endo and exo conformational products were shown. The second-order reaction rate constant k2 is 4.4 × 10−2 M−1 s−1 for endo pathway and 0.18 × 10−2 M−1 s−1 for exo pathway, with a ratio of 25 from the experimental measurement.

exo-cis conformation.17 This −0.44 kcal mol−1 difference in the reaction barrier still deviates from the experimental value (−1.9 kcal mol−1). Moreover, the absolute reaction barriers for both reaction pathways deviate significantly from the experimental ones, leading to incorrect absolute reaction rates for this reaction. To a certain extent, there was a large error cancellation in the calculation of endo/exo selectivity for this DA reaction. In their work,17 both CP and MVK were optimized at RHF/6-31G* level in the gas phase. The energies were calculated at MP3/6-31G*//6-31G* level, and the vibration frequencies for all reactants and transition states were scaled by a factor of 0.91. They found that the calculated entropies had a constant deviation of about 8 kcal mol−1 K−1 with those obtained from the experiment, which might be corrected by an expected translational entropy in solution. Later on, Jorgensen and co-workers studied the solvation effect on this reaction by using the AM1/MM/Monte Carlo method with a one-dimensional (1D) reaction coordinate. They found that the hydration effect on the reaction rate was consistent with experimental observations and the observed increase in the reaction rate in water was primarily originated from enhanced hydrogen bonding between the solvent and the polar transition state.19 They also noticed that a single reaction coordinate fails to provide an accurate location of transition structures. To get around this problem, Jorgensen and coworkers continued to study this reaction by the PDDG/PM3/ MM/Monte Carlo method with two-dimensional (2D) reaction coordinates, where the acceleration of this DA reaction in water is further verified by the enhanced hydrogen bonding to hydrogen-bond-accepting groups in dienophile.20,24 However, theoretical results overestimated the experimental reaction barrier by about 10 kcal mol−1, unless corrections at MP2/6-311+G(d,p) level were applied. Recently, Houk and co-workers developed a solvent-perturbed transition-state (SPTS) sampling scheme to simulate chemical reaction dynamics in the condensed phase and applied this method to this DA reaction in water.21 Their results supported the observations of Jorgensen and co-workers and indicated hydrogen-bond formation between the solvent and the TS by reorientation of a solvent molecule during the reaction process. On the other hand, the authors suggested that the polarization from the water molecules surrounding the TS is essential for the observation of hydrogen-bond reorientation of water molecules adjacent to the reacting molecules.22 However, an accurate calculation of the reaction barrier was not considered in their study. To give a reliable explanation to the stereoselectivity of this DA reaction, a free-energy surface at an ab initio (ai) level is indispensable. However, the direct computation of a free-

energy surface at an ai level is formidably expensive, even if it is not impossible. Recently, one variant of the reference-potential methods25,26 termed MBAR+wTP27 was proposed by some of the authors, via which an accurate ai quantum mechanical/ molecular mechanical (QM/MM) free-energy (FE) profile can be computed with much less computational expense. In this method, a weighted thermodynamic perturbation (wTP)28 correction is applied to an inexpensive profile, which is generated by the multistate Bennett acceptance ratio (MBAR)29,30 analysis of the trajectories from umbrella samplings (US).31 This MBAR+wTP method has been successfully applied to calculate the free-energy profiles of one quasichemical reaction and three chemical reactions in aqueous solution.27 We also applied this method to the calculations of the free-energy surfaces of two Diels−Alder reactions of cyclopentadiene with either acrylonitrile or 1-4naphthoquinone.32 In the current work, we will investigate the origin of this endo/exo selectivity of the DA reaction between CP and MVK in aqueous solution at ai QM/MM level in an explicit solvent using the MBAR+wTP method. Comparison with routinely applied stationary structure methods in the implicit solvent will also be provided.

■

METHODS Multistate Bennett Acceptance Ratio and Weighted Thermodynamic Perturbation (MBAR+wTP) Method. The MBAR+wTP method has been explained at length in ref 27. Here, we will only briefly review its main idea. In the first step, umbrella sampling (US) simulations under an inexpensive Hamiltonian, e.g., semiempirical (SE) QM/MM, are performed. Overlap in the phase space between two neighboring windows is essential for a subsequent analysis. Therefore, it is monitored by the overlap matrix proposed by Klimovich et al.;33 then, the FE profile at this level of theory is calculated using the multistate Bennett acceptance ratio (MBAR) analysis; finally, weighted thermodynamic perturbation (wTP) from the SE Hamiltonian to the ai Hamiltonian is performed to obtain the ai QM/MM FE profile using weight factors from the previous MBAR analysis. The reliability of this wTP calculation is characterized by the reweighting entropy,34 and Gaussian process regression (GPR) method35 is utilized to remove the statistical noise in the wTP correction and smooth the FE profile. Structural-ensemble averages of any physical properties can be computed with these weight factors and the properties for each sample. Calculation of a 1D FE Profile from a 2D FE Surface. After a two-dimensional (2D) free-energy (FE) surface is obtained by the MBAR+wTP method, the curvilinear 5132

DOI: 10.1021/acs.jpcb.9b01989 J. Phys. Chem. B 2019, 123, 5131−5138

Article

The Journal of Physical Chemistry B

Figure 2. Free-energy surfaces under the PM6 Hamiltonian (left) and under B3LYP Hamiltonian after Gaussian process regression (right) for (a) the endo pathway and (b) the exo pathway. The minimum free-energy paths are shown as red-dotted lines.

where ΔG‡ is the activation free energy, kB is the Boltzmann constant, T is the absolute temperature, and h is the Planck constant. κ is the transmission coefficient, which is assumed to be equal to one, assuming no recrossing takes place. Thus, the ratio of the endo and exo reaction rate constants can be calculated from the difference in the activation free energies of these two pathways as É ÅÄÅ ‡ ‡ Ñ ÑÑ ÅÅ ΔGendo − ΔGexo kendo ÑÑ Å = expÅÅ− rendo/exo = ÑÑ Å ÑÑ kexo kBT ÅÅÇ (5) ÑÖ

minimum free-energy path (MFEP) on the free-energy surface can be located by using MEPSA software.36 Along the MFEP, a 1D FE profile can be calculated by the 2D FE profile as follows. After grouping all grids on the 2D FE surface with free energies f i(η1, η2) to its nearest point with a generalized reaction coordinate s on the MFEP, the 1D FE profile can be calculated by G(s) = −kBT ln P(s)

in which

i

fi yz zz zz k T B k {

∑ expjjjjj− Ns

P(s) =

(1)

i=1

Computational Details. As shown in Figure 1, the Diels− Alder reaction between cyclopentadiene (CP) and methyl vinyl ketone (MVK) was studied in this work, where s-cis conformation for MVK was used and both the endo and exo addition modes for this reaction are displayed. The s-cis conformation for MVK and endo addition mode is preferred from experimental investigations12,16 and is further evidenced by some theoretical calculations on the energetics of the transition structures.17,20 The second-order reaction rate k2 for endo product is 4.4 × 10−2 M−1 s−1, and the experimentally measured ratio of the endo/exo products is about 25.10,12,13,15 In this work, both CP and MVK were included in the SQM or QM region, which was surrounded by a TIP3P water sphere with a radius of 25 Å centering on the heavy atom closest to the center of mass of the SQM or QM region. Although simulations under periodic boundary condition are more rigorous than that in a water droplet,38−42 Vasilevskaya et al. have verified that the periodic and nonperiodic QM/MM treatments can give similar free-energy profiles for the reactions in solution.43 Their work also shows that the longrange electrostatic interactions can be well captured by nonperiodic QM/MM calculations in a water droplet of

(2)

Here, Ns is the total number of all points on the 2D FE surface belonging to the 1D reaction coordinate s on the MFEP according to the principle of nearest distance redistribution. From the error propagation rule, the variance for the 1D FE profile can be computed via N

σG2

2

= (kBT )

∑i =s 1

2

( )σ xi kBT

N

(∑i =s 1 xi)2

2 fi

(3)

f

( )

where xi = exp − k Ti and σ2f i is the variance of each grid on B

the FE surface and can be obtained by GPR, as mentioned above. With the 1D FE profile, the reaction rate constant k2 can be computed according to the Eyring equation37 ÄÅ ÉÑ κkBT ÅÅÅ ΔG‡ ÑÑÑ Å ÑÑ k2 = expÅÅ− ÅÅÇ kBT ÑÑÑÖ h (4) 5133

DOI: 10.1021/acs.jpcb.9b01989 J. Phys. Chem. B 2019, 123, 5131−5138

Article

The Journal of Physical Chemistry B reasonable size.43 Therefore, the droplet model is widely used in QM/MM studies of enzymatic reactions and chemical reactions in solution. This water sphere was restrained by a soft half-harmonic potential with a force constant of 10 kcal mol−1 Å−2 to avoid evaporation. The nonbonded interactions were fully counted without any truncation. The van der Waals (vdW) parameters for the SQM and QM regions were taken from the general AMBER force field.44 PM6 and B3LYP with a basis set of 6-31G* were chosen as the low-level and the highlevel Hamiltonians, respectively. Two-dimensional umbrella sampling (US) simulations were performed at PM6 level, and the indirect FE surface at B3LYP/MM level was obtained by using our MBAR+wTP method. For the stationary states (reactant, product, and transition states), the geometry optimization and frequency analysis were also performed at B3LYP and the 2nd order Møller−Plesset perturbation theory (MP2) levels using a basis set of 6-31G* and the polarizable continuum solvent model using the integral equation formalism variant (IEFPCM) for the solvation effect. The two-dimensional US windows center on η ≡ (η1 = dC1C3, η2 = dC2C4) ranging from 1.50 to 4.00 Å with an increment of 0.05 Å for each dimension. The reactant state was defined as η1 = 4.00 Å and η2 = 4.00 Å, and the product state was defined as η1 = 1.58 Å and η2 = 1.58 Å. In each US window simulation under the PM6/MM Hamiltonian, after a minimization in energy, the system was heated up to 298.15 K in 50 ps and then was equilibrated for 100 ps. At last, a 2 ns production simulation was performed for a subsequent freeenergy analysis, in which the integration time step was set to 1 fs and configurations were saved every 1 ps. The Andersen-like temperature coupling scheme45 was used to regulate the temperature at 298.15 K. Then, single-point energies under the PM6/MM and B3LYP/MM Hamiltonians for 1000 configurations were computed for the wTP calculations. MEPSA software36 was used to search for the minimum free-energy paths on the free-energy surfaces. AmberTools 17 program package46 was used to perform all US simulations and the single-point energy calculations, by interfacing with Gaussian 16 package47 when the energies at B3LYP level were requested. The calculations in an implicit solvent were also carried out using Gaussian 16 package.

Figure 3. Free-energy profiles for the endo (red) and exo (green) reactions along the minimum free-energy path under PM6 Hamiltonian (top) and under B3LYP Hamiltonian (bottom). Insets are the free-energy profiles near the transition states.

As listed in Table 1, Jorgensen et al. showed that the activation free energy and reaction free energy for this reaction are 32.2 ± 0.5 and −25.8 ± 0.6 kcal mol−1, respectively, under

■

RESULTS DA Reaction toward an Endo Isomer between CP and MVK. As shown in Figure 2, for this DA reaction with the endo isomer, the structure with dC1C3 = 4.00 Å and dC2C4 = 4.00 Å is defined as the reactant and that with dC1C3 = 1.58 Å and dC2C4 = 1.58 Å is defined as the product. The transition state can be located at (η1, η2) = (2.47, 1.88 Å) at PM6 level, which has some deviation from (η1, η2) = (2.34, 2.01 Å) at PM3 level.20 The transition state at B3LYP level is located at (η1, η2) = (2.57, 2.10 Å), showing longer bond lengths than those in the transition state at PM6 level but still indicating apparent asynchronism of this reaction. The transition states at both the PM6 and B3LYP levels are located at the position of the maximum along their respective minimum free-energy paths, as shown in Figure 3. In the implicit solvent model (PCM), this asynchronism is further manifested at the B3LYP level, whose transition state is located at (η1, η2) = (2.66, 2.02 Å). The MP2 method in the implicit solvent model locates the transition state at (η1, η2) = (2.52, 2.19 Å), showing that small deviations from the B3LYP/MM level result by less than 0.1 Å.

Table 1. Activation Free Energy ΔG‡ and Reaction Free Energy ΔG (in a Unit of kcal mol−1) at 298.15 K product

method

endo

PM3b AM1c PM6 B3LYP B3LYP/PCM MP2/PCM exp.d PM6 B3LYP B3LYP/PCM MP2/PCM exp.d

exo

ΔG‡ 32.2 22.6 30.7 18.6 28.9 14.2 19.3 30.0 19.1 29.4 15.8 21.2

± 0.5 ± 0.1 ± 0.3

± 0.2 ± 0.5

ΔG

(η1, η2)a

−25.8 ± 0.6 −24.6 −16.0 ± 0.1 −13.1 ± 0.2 −0.07 −14.2

(2.34, 2.01) (2.47, (2.57, (2.66, (2.52,

1.88) 2.10) 2.02) 2.19)

−16.6 ± 0.1 −13.2 ± 0.2 −0.3 −13.6

(2.45, (2.56, (2.62, (2.48,

1.91) 2.05) 2.03) 2.20)

a

Transition-state position in Å. bRef 20 by Jorgensen et al. using PDDG/PM3/MM/MC. cRef 19 with one-dimensional reaction coordinate. dRef 15.

5134

DOI: 10.1021/acs.jpcb.9b01989 J. Phys. Chem. B 2019, 123, 5131−5138

Article

The Journal of Physical Chemistry B PDDG/PM3 Hamiltonian,20 and 22.6 and −24.6 kcal mol−1 under AM1 Hamiltonian.19 The 1D FE profiles at PM6 and B3LYP levels for this DA reaction are shown in Figure 3 (red curve). The results show a strong dependence on Hamiltonians. The reaction barrier and the reaction free energy are 30.7 ± 0.1 and −16.0 ± 0.1 kcal mol−1, respectively, under PM6 Hamiltonian. The activation free energy of this reaction, which has an experimental value of 19.3 kcal mol−1,15 was significantly overestimated at the PM6 level. This large overestimate of the activation barrier originates from the approximations lying in these semiempirical methods20 and is also observed from two similar reactions in our previous work.32 An activation barrier of 18.6 ± 0.3 kcal mol−1 at the B3LYP/MM level is obtained by our MBAR+wTP method, which is much closer to the experimental value than those obtained by the semiempirical methods, including the PM3 and AM1 potentials. The reliability of the wTP calculation can be characterized by the reweighting entropy, as mentioned in our previous works.27,32 This indicates that the phase space overlap between the PM6/MM and B3LYP/MM Hamiltonians is sufficient. Therefore, the wTP calculations and the weighted ensemble-average calculations are reliable. It should be noted that the B3LYP/PCM method overestimates the activation free energy by 9.6 kcal mol−1 whereas the MP2/ PCM method underestimates the activation free energy by 5.1 kcal mol−1, which also agrees with what has been observed in our previous work.32 Shown in Figure 4 are the radial

PCM level, compared with −24.6 kcal mol−1 at AM1 level, −25.8 ± 0.6 kcal mol−1 at PM3 level, and −16.0 ± 0.1 kcal mol−1 at PM6 level. However, B3LYP/PCM remarkably underestimates the exothermicity of this reaction and shows almost no exothermicity for this reaction. DA Reaction toward an Exo Isomer between CP and MVK. From the 2D free-energy landscape shown in Figure 2, the reactant and the product of this exo pathway are defined similarly as in the endo pathway shown above. The transition state can be located at (η1, η2) = (2.45, 1.91 Å) at PM6 level, with η1 slightly shorter but η2 slightly longer than their counterparts in the endo pathway. Again, the B3LYP/MM potential prefers longer distances for the reaction coordinates ((η1, η2) = (2.56, 2.05 Å)) at the transition states than the PM6 potential and it also exhibits the asynchronism of this reaction. At the B3LYP/PCM level, the transition state is located at (η1, η2) = (2.62, 2.03 Å). With MP2 potential replacing B3LYP, the transition state in the implicit solvent model becomes shorter in η1 by 0.14 Å and longer in η2 by 0.17 Å. Overall, the transition structures of these two reaction pathways are quite similar regardless of the Hamiltonians, at least for those used in this study, with the difference in bond lengths on a scale of 0.01 Å. The 1D FE profiles at PM6 and B3LYP levels for this DA reaction are shown in Figure 3 (green curve). The activation free energy and reaction free energy for this reaction are 30.0 ± 0.2 and −16.6 ± 0.1 kcal mol−1, respectively, under PM6 Hamiltonian, differing by only 0.6−0.7 kcal mol−1 from those of the endo reaction pathway under the same Hamiltonian, as shown in Table 1. The almost identical reaction barriers can be verified by the high similarity in the RDF curves, as shown in the top panel of Figure 4. This semiempirical method still overestimates the activation barrier of this reaction, which has an experimental value of 21.2 kcal mol−1.15 Our MBAR+wTP method yields an activation barrier of 19.1 ± 0.5 kcal mol−1 at the B3LYP/MM level, which is much closer to the experimental value than that at PM6/MM level. The activation free energy is overestimated by 8.2 kcal mol−1 at B3LYP/PCM level and is underestimated by 5.4 kcal mol−1 at MP2/PCM level. These are very similar to those in the endo pathway. As shown in the bottom panel of Figure 4, under B3LYP Hamiltonian, the mean RDF curve between the oxygen atom in MVK and the hydrogen atoms in water molecules in the TS ensemble of this exo pathway has a smaller g(R) value and a slightly longer distance at the first peak than those for the endo pathway, which can explain the higher reaction barrier of the exo pathway than that of the endo one. Meanwhile, the B3LYP/MM method gives a reaction free-energy value of −13.2 ± 0.2 kcal mol−1, which agrees much better with MP2/ PCM (−13.6 kcal mol−1) than the PM6 potential (−16.6 ± 0.1 kcal mol−1). B3LYP/PCM still remarkably underestimates the exothermicity of this reaction. As shown in Figure 5, we chose the one with the maximum MBAR weight from the many configurations in the bin holding the transition-state geometry at the B3LYP/MM level. Both configurations could be approximately treated as the transitionstate structures at B3LYP/MM level for the endo pathway and the exo pathway, respectively. For the endo pathway, there are three water molecules forming hydrogen bonds with the oxygen atom in MVK whereas for the exo pathway, there are only two water molecules hydrogen-bonded to MVK. This agrees well with the observation that the mean RDF curve between the oxygen atom in MVK and the hydrogen atoms in

Figure 4. Solvent radial distribution near the oxygen atom in MVK averaged over the transition-state ensembles under PM6 Hamiltonian (top) and under B3LYP Hamiltonian (bottom) for the endo pathway (red) and the exo pathway (green).

distribution functions (RDF) averaged over the TS ensembles for this endo pathway under PM6 and B3LYP levels of theory (red curve). The larger g(R) value and the shorter distance between the oxygen atom in MVK and the hydrogen atom in the solvent at the first peak of the RDF curve under the B3LYP Hamiltonian indicate stronger hydrogen bonding interaction between water molecules and the TS structure, compared to those under the PM6 Hamiltonian, which explains the lower reaction barrier at the B3LYP level and is also consistent with the observations of Jorgensen and co-workers.20 Meanwhile, the B3LYP/MM level used in our MBAR+wTP method can give a reaction free-energy value of −13.1 ± 0.2 kcal mol−1, which is much closer to the value of −14.2 kcal mol−1 at MP2/ 5135

DOI: 10.1021/acs.jpcb.9b01989 J. Phys. Chem. B 2019, 123, 5131−5138

Article

The Journal of Physical Chemistry B

Figure 5. Near-transition-state structures at B3LYP/MM level for the endo pathway (left) and exo pathway (right), respectively. The structures were taken from the bins holding the transition state with the maximum weight at the B3LYP/MM level.

water molecules in the TS ensemble of this exo pathway has a smaller g(R) value and a slightly longer distance at the first peak than those for the endo pathway, which can explain the higher reaction barrier of the exo pathway than that of the endo one. Endo/Exo Selectivity of the DA Reaction between CP and MVK. As shown in Table 2, the experimental reaction rate

Table 3. Estimated Wall-Clock Time in Hours for Computing the QM/MM Free-Energy Surfaces at B3LYP/631G* Levela direct B3LYP/MMb

MBAR+wTP sampling

energy evaluation

total

total

5800

1800

7600

1 500 000

a

Assuming one node with 16 cores of Intel Xeon CPU E5-2660 2.20 GHz was used. bUsing the same number of windows as in the semiempirical simulations, and one 1 ns simulation for each window.

Table 2. Reaction Rate Constants (in the Unit M−1 s−1) for the Endo and Exo Pathways as Well as the Ratio of Endo and Exo Products method PM6 B3LYP B3LYP/PCM MP2/PCM exp.a

kendo 2.0 × 1.4 × 4.1 × 242.5 4.4 ×

kexo −10

10 10−1 10−9 10−2

6.4 × 6.2 × 1.8 × 16.3 1.8 ×

landscape at the B3LYP/MM level. Therefore, the MBAR +wTP method is about 200 times faster than a direct simulation.

rendo/exo −10

10 10−2 10−9 10−3

0.3 2 2 15 25

■

DISCUSSION For studies of the reaction mechanism and thus calculations of the free-energy landscape at a high level of theory, the method of choice must maintain a good balance between the sampling efficiency and statistical efficiency. Although a direct calculation of the free-energy landscape at the ai level has high statistical efficiency, it has very poor sampling efficiency. The typical simulation time is limited to tens of picoseconds. On the contrary, reference-potential methods have high sampling efficiency but they suffer from low statistical efficiency due to exponential terms in the low-level to highlevel correction. This statistical inefficiency is originated from the insufficient overlap in the phase space between these two Hamiltonians.34,48−51 To increase the overlap, we can either modify the low-level Hamiltonian (e.g., via force matching52) or supplement a correction term (e.g., represented as an artificial neural network53) to the low-level Hamiltonian. In this work, this small difference in the free-energy barrier (1.9 kcal mol−1) poses a challenge to the computational methods, which can be easily submerged in statistical uncertainties. In this regard, a more robust free-energy calculation strategy or pathway is indispensable. In this study, the vdW interaction between the (S)QM region and the MM region is computed using molecular mechanics. Besides the empirical spirit of these vdW parameters, the atom types, using the language of force fields, are not fixed during the reaction. For instance, four sp2 carbon atoms become sp3 atoms in this DA reaction. Their vdW parameters should also be changed accordingly. A more “ab initio” description of the dispersion term should implement

a

Refs 10, 12, 13, 15.

constants for the endo and exo pathways are 4.4 × 10−2 and 1.8 × 10−3 M−1 s−1, respectively, which leads to an endo/exo ratio of about 25. Due to the approximations in the PM6 Hamiltonian, both of the calculated reaction rate constants deviate remarkably from the experimental measurements, which results in an endo/exo ratio of 0.3 with a preference to the exo product. However, after reweighting from the PM6/ MM results, the B3LYP/MM method yields two reaction freeenergy barriers that are much closer to the experimental values than those obtained at PM6, B3LYP/PCM, and MP2/PCM levels. The endo/exo ratio from the B3LYP/MM method is 2, which is better than the one at the PM6/MM level but is still 1 order smaller than the experimental one. Although MP2/PCM shows the best agreement with the experiment in terms of the endo/exo ratio (15), the absolute reaction free-energy barriers show deviations of about 5 kcal mol−1 and the success in the prediction of the stereoselectivity is a fortuitous error cancellation. CPU Time. The CPU time for the calculations is listed in Table 3. Using the MBAR+wTP method, the dominant CPU time is spent on the umbrella sampling simulations, which is 5800 h, whereas the single-point calculations at PM6/MM and B3LYP/MM levels take only 1800 h. In total, the CPU time is about 7600 h for each reaction pathway. However, it takes about 1 500 000 h for direct calculation of the free-energy 5136

DOI: 10.1021/acs.jpcb.9b01989 J. Phys. Chem. B 2019, 123, 5131−5138

Article

The Journal of Physical Chemistry B dynamic atomic dipole polarizabilities during the reaction, which is a function of the electron density partitioned to the specific atom theoretically.54−56 However, this may significantly increase the complexity of the model, and some other difficulty may also arise, for instance, the nonunique partition of the electron density onto each atom.57 In this study, the reaction coordinates were chosen intuitively as the lengths of the newly formed C−C bonds. On the one hand, the mapping of a two-dimensional freeenergy landscape to a one-dimensional one is not well defined, especially as in this case, the minimum free-energy path in the two-dimensional free-energy landscape is a nonlinear curve. The directions perpendicular to the tangents at two points on the minimum free-energy path may cross somewhere off the path, which leads to ambiguity in integration regions when doing the 2D-to-1D mapping. In this work, we grouped the grids by their distances to the discrete points on the minimum free-energy path, which is actually a Voronoi diagram. However, this is never the unique way to do this 2D-to-1D mapping. On the other hand, the selection of the reaction coordinate is nontrivial and is still an open topic. An optimal choice of the reaction coordinate should be able to represent the committor hypersurface as well as be possible especially for the transition state but should be as low as possible in dimension (least redundant). However, evaluation of the quality of the reaction coordinate may lead to an extremely high computational cost at the ai QM/MM level.58

Ye Mei: 0000-0002-3953-8508 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Y.M. is supported by the Ministry of Science and Technology of China (Grant No. 2016YFA0501700), the National Natural Science Foundation of China (Grant No. 21773066), and the Fundamental Research Funds for the Central Universities. Y.S. is grateful to Department of Energy Office of Science (DESC0011297). CPU time was provided by the Supercomputer Center of East China Normal University (ECNU Multifunctional Platform for Innovation No. 001).

■

(1) Kloetzel, M. C. Organic Reactions; American Cancer Society, 2011; Chapter 1, pp 1−59. (2) Holmes, H. L. Organic Reactions; American Cancer Society, 2011; Chapter 2, pp 60−173. (3) Carruthers, W. Cycloaddition Reactions in Organic Synthesis; Tetrahedron Organic Chemistry Series; Elsevier, 1990; Vol. 8, pp 91− 139. (4) Nicolaou, K. C.; Snyder, S. A.; Montagnon, T.; Vassilikogiannakis, G. The Diels-Alder Reaction in Total Synthesis. Angew. Chem., Int. Ed. 2006, 41, 1668−1698. (5) Loncharich, R. J.; Brown, F. K.; Houk, K. N. Transition Structures of the Diels-Alder Reaction of Butadiene with Acrolein. J. Org. Chem. 1989, 54, 1129−1134. (6) Birney, D. M.; Houk, K. N. Transition Structures of the Lewis Acid-Catalyzed Diels-Alder Reaction of Butadiene with Acrolein. The Origins of Selectivity. J. Am. Chem. Soc. 1990, 112, 4127−4133. (7) Bernardi, F.; Bottoni, A.; Field, M. J.; Guest, M. F.; Hillier, I. H.; Robb, M. A.; Venturini, A. J. MC-SCF Study of the Diels-Alder Reaction between Ethylene and Butadiene. J. Am. Chem. Soc. 1988, 110, 3050−3055. (8) Li, Y.; Houk, K. N. Diels-Alder Dimerization of 1,3-Butadiene: An ab initio CASSCF Study of the Concerted and Stepwise Mechanisms and Butadiene-Ethylene Revisited. J. Am. Chem. Soc. 1993, 115, 7478−7485. (9) Sustmann, R.; Sicking, W. Influence of Reactant Polarity on the Course of (4 + 2) Cycloadditions. J. Am. Chem. Soc. 1996, 118, 12562−12571. (10) Rideout, D. C.; Breslow, R. Hydrophobic Acceleration of DielsAlder Reactions. J. Am. Chem. Soc. 1980, 102, 7816−7817. (11) Breslow, R.; Maitra, U.; Rideout, D. C. Selective Diels-Alder Reactions in Aqueous Solutions and Suspensions. Tetrahedron Lett. 1983, 24, 1901−1904. (12) Breslow, R.; Maitra, U. On the Origin of Product Selectivity in Aqueous Diels-Alder Reactions. Tetrahedron Lett. 1984, 25, 1239− 1240. (13) Breslow, R.; Guo, T. Diels-Alder Reactions in Nonaqueous Polar Solvents. Kinetic Effects of Chaotropic and Antichaotropic Agents and of β-cyclodextrin. J. Am. Chem. Soc. 1988, 110, 5613− 5617. (14) Breslow, R. Hydrophobic Effects on Simple Organic Reactions in Water. Acc. Chem. Res. 1991, 24, 159−164. (15) Engberts, J. B. Diels-Alder Reactions in Water: Enforced Hydrophobic Interaction and Hydrogen Bonding. Pure Appl. Chem. 1995, 67, 823−828. (16) Cativiela, C.; Garcia, J.; Gil, J.; Martinez, R.; Mayoral, J.; Salvatella, L.; Urieta, J.; Mainar, A.; Abraham, M. Solvent Effects on Diels-Alder Reactions. The Use of Aqueous Mixtures of Fluorinated Alcohols and the Study of Reactions of Acrylonitrile. J. Chem. Soc., Perkin Trans. 2 1997, 3, 653−660. (17) Jorgensen, W. L.; Lim, D.; Blake, J. F. Ab initio Study of DielsAlder Reactions of Cyclopentadiene with Ethylene, Isoprene,

■

CONCLUSIONS In this work, the endo/exo stereoselectivity of the Diels−Alder (DA) reaction between cyclopentadiene (CP) and methyl vinyl ketone (MVK) in an aqueous solution has been investigated using a method termed MBAR+wTP proposed by us previously. This method is based on the main idea of the reference-potential methods. Accurate absolute free-energy barriers at B3LYP/6-31G* level for these two reactions in explicit water molecules were computed indirectly with muchenhanced efficiency, which are very close to the experimental measurements. The stereoselectivity mainly comes from the solvation effect to the polar transition state. The first peak of the solvent distribution near the oxygen atom in MVK in the transition-state ensemble is slightly closer for the endo pathway than that in the exo pathway. Nonetheless, we obtained an endo/exo ratio of 2 from the difference of the reaction freeenergy barrier at the B3LYP/MM level, which is still 1 order of magnitude smaller than the experimental measurement and the uncertainty is too large. To obtain more reliable stereoselectivity of this reaction, a more robust free-energy method should be employed, for instance, a direct conversion between two transition ensembles. Relatively speaking, the semiempirical method PM6/MM significantly overestimates the reaction free-energy barriers. Although the endo/exo ratio from the MP2/PCM calculation agrees well with the experiment, the absolute reaction barriers are off by about 5 kcal mol−1.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Pengfei Li: 0000-0001-9355-6371 Yihan Shao: 0000-0001-9337-341X 5137

DOI: 10.1021/acs.jpcb.9b01989 J. Phys. Chem. B 2019, 123, 5131−5138

Article

The Journal of Physical Chemistry B Cyclopentadiene, Acrylonitrile, and Methyl Vinyl Ketone. J. Am. Chem. Soc. 1993, 115, 2936−2942. (18) Blake, J. F.; Lim, D.; Jorgensen, W. L. Enhanced Hydrogen Bonding of Water to Diels-Alder Transition States. ab initio Evidence. J. Org. Chem. 1994, 59, 803−805. (19) Chandrasekhar, J.; Shariffskul, S.; Jorgensen, W. L. QM/MM Simulations for Diels-Alder Reactions in Water: Contribution of Enhanced Hydrogen Bonding at the Transition State to the Solvent Effect. J. Phys. Chem. B 2002, 106, 8078−8085. (20) Acevedo, O.; Jorgensen, W. L. Understanding Rate Accelerations for Diels-Alder Reactions in Solution Using Enhanced QM/MM Methodology. J. Chem. Theory Comput. 2007, 3, 1412− 1419. (21) Yang, Z.; Doubleday, C.; Houk, K. N. QM/MM Protocol for Direct Molecular Dynamics of Chemical Reactions in Solution: The Water-Accelerated Diels-Alder Reaction. J. Chem. Theory Comput. 2015, 11, 5606−5612. (22) Liu, F.; Yang, Z.; Mei, Y.; Houk, K. N. QM/QM′ Direct Molecular Dynamics of Water-Accelerated Diels-Alder Reaction. J. Phys. Chem. B 2016, 120, 6250−6254. (23) Soto-Delgado, J.; Tapia, R. A.; Torras, J. Multiscale Treatment for the Molecular Mechanism of a Diels-Alder Reaction in Solution: A QM/MM-MD Study. J. Chem. Theory Comput. 2016, 12, 4735−4742. (24) Thomas, L. L.; Tirado-Rives, J.; Jorgensen, W. L. Quantum Mechanical/Molecular Mechanical Modeling Finds Diels-Alder Reactions Are Accelerated Less on the Surface of Water Than in Water. J. Am. Chem. Soc. 2010, 132, 3097−3104. (25) Gao, J. Absolute Free Energy of Solvation from Monte Carlo Simulations Using Combined Quantum and Molecular Mechanical Potentials. J. Phys. Chem. 1992, 96, 537−540. (26) Muller, R. P.; Warshel, A. Ab Initio Calculations of Free Energy Barriers for Chemical Reactions in Solution. J. Phys. Chem. 1995, 99, 17516−17524. (27) Li, P.; Jia, X.; Pan, X.; Shao, Y.; Mei, Y. Accelerated Computation of Free Energy Profile at ab initio QM/MM Accuracy via a Semi-Empirical Reference-Potential. I. Weighted Thermodynamics Perturbation. J. Chem. Theory Comput. 2018, 14, 5583. (28) Zwanzig, R. W. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases. J. Chem. Phys. 1954, 22, 1420. (29) Shirts, M. R.; Chodera, J. D. Statistically Optimal Analysis of Samples from Multiple Equilibrium States. J. Chem. Phys. 2008, 129, No. 124105. (30) Shirts, M. R. Reweighting from the Mixture Distribution as a Better Way to Describe the Multistate Bennett Acceptance Ratio. 2017, arXiv:1704.00891. arXiv.org e-Print archive. https://arxiv.org/ abs/1704.00891. (31) Torrie, G.; Valleau, J. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23, 187−199. (32) Li, P.; Liu, F.; Jia, X.; Shao, Y.; Hu, W.; Zheng, J.; Mei, Y. Efficient Computation of Free Energy Surfaces of Diels-Alder Reactions in Explicit Solvent at ab initio QM/MM Level. Molecules 2018, 23, No. 2487. (33) Klimovich, P. V.; Shirts, M. R.; Mobley, D. L. Guidelines for the Analysis of Free Energy Calculations. J. Comput.-Aided Mol. Des. 2015, 29, 397−411. (34) Wang, M.; Li, P.; Jia, X.; Liu, W.; Shao, Y.; Hu, W.; Zheng, J.; Brooks, B. R.; Mei, Y. Efficient Strategy for the Calculation of Solvation Free Energies in Water and Chloroform at the Quantum Mechanical/Molecular Mechanical Level. J. Chem. Inf. Model. 2017, 57, 2476−2489. (35) Rasmussen, C.; Williams, C. Gaussian Processes for Machine Learning; MIT Press, 2006; pp 7−32. (36) Marcos-Alcalde, I.; Setoain, J.; Mendieta-Moreno, J. I.; Mendieta, J.; Gómez-Puertas, P. MEPSA: Minimum Energy Pathway Analysis for Energy Landscapes. Bioinformatics 2015, 31, 3853−3855. (37) Wynne-Jones, W. F. K.; Eyring, H. The Absolute Rate of Reactions in Condensed Phases. J. Chem. Phys. 1935, 3, 492−502.

(38) Nam, K.; Gao, J.; York, D. M. An Efficient Linear-Scaling Ewald Method for Long-Range Electrostatic Interactions in Combined QM/ MM Calculations. J. Chem. Theory Comput. 2005, 1, 2−13. (39) Walker, R. C.; Crowley, M. F.; Case, D. A. The Implementation of a Fast and Accurate QM/MM Potential Method in Amber. J. Comput. Chem. 2008, 29, 1019−1031. (40) Giese, T. J.; York, D. M. Ambient-Potential Composite Ewald Method for ab Initio Quantum Mechanical/Molecular Mechanical Molecular Dynamics Simulation. J. Chem. Theory Comput. 2016, 12, 2611−2632. (41) Holden, Z. C.; Rana, B.; Herbert, J. M. Analytic Gradient for the QM/MM-Ewald Method using Charges Derived from the Electrostatic Potential: Theory, Implementation, and Application to ab initio Molecular Dynamics Simulation of the Aqueous Electron. J. Chem. Phys. 2019, 150, No. 144115. (42) Kawashima, Y.; Ishimura, K.; Shiga, M. Ab initio Quantum Mechanics/Molecular Mechanics Method with Periodic Boundaries Employing Ewald Summation Technique to Electron-Charge Interaction: Treatment of the Surface-Dipole Term. J. Chem. Phys. 2019, 150, No. 124103. (43) Vasilevskaya, T.; Thiel, W. Periodic Boundary Conditions in QM/MM Calculations: Implementation and Tests. J. Chem. Theory Comput. 2016, 12, 3561. (44) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (45) Andersen, H. C. Molecular Dynamics Simulations at Constant Pressure and/or Temperature. J. Chem. Phys. 1980, 72, 2384−2393. (46) Case, D. A.; Cerutti, D. S.; Cheatham, T. E.; Darden, I. I. I.; Duke, T. A.; Giese, R. E.; Gohlke, T. J.; Goetz, H.; Greene, A. W.; Homeyer, D. AMBER 2017; University of California: San Francisco, 2017. (47) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H. Gaussian 16; Gaussian Inc.: Wallingford, CT, 2016. (48) Kofke, D. A. On the Sampling Requirements for ExponentialWork Free-Energy Calculations. Mol. Phys. 2006, 104, 3701−3708. (49) Pohorille, A.; Jarzynski, C.; Chipot, C. Good Practices in FreeEnergy Calculations. J. Phys. Chem. B 2010, 114, 10235−10253. (50) Boresch, S.; Woodcock, H. L. Convergence of Single-Step Free Energy Perturbation. Mol. Phys. 2017, 115, 1200−1213. (51) Ryde, U. How Many Conformations Need to Be Sampled to Obtain Converged QM/MM Energies? The Curse of Exponential Averaging. J. Chem. Theory Comput. 2017, 13, 5745−5752. (52) Zhou, Y.; Pu, J. Reaction Path Force Matching: A New Strategy of Fitting Specific Reaction Parameters for Semiempirical Methods in Combined QM/MM Simulations. J. Chem. Theory Comput. 2014, 10, 3038−3054. (53) Shen, L.; Yang, W. Molecular Dynamics Simulations with Quantum Mechanics/Molecular Mechanics and Adaptive Neural Networks. J. Chem. Theory Comput. 2018, 14, 1442−1455. (54) Tkatchenko, A.; DiStasio, R. A.; Car, R.; Scheffler, M. Accurate and Efficient Method for Many-Body van der Waals Interactions. Phys. Rev. Lett. 2012, 108, No. 236402. (55) Kuechler, E. R.; Giese, T. J.; York, D. M. Charge-Dependent Many-Body Exchange and Dispersion Interactions in Combined QM/ MM Simulations. J. Chem. Phys. 2015, 143, No. 234111. (56) Caldeweyher, E.; Ehlert, S.; Hansen, A.; Neugebauer, H.; Spicher, S.; Bannwarth, C.; Grimme, S. A Generally Applicable Atomic-Charge Dependent London Dispersion Correction. J. Chem. Phys. 2019, 150, No. 154122. (57) Mei, Y.; Simmonett, A. C.; Pickard, F. C.; DiStasio, R. A.; Brooks, B. R.; Shao, Y. Numerical Study on the Partitioning of the Molecular Polarizability into Fluctuating Charge and Induced Atomic Dipole Contributions. J. Phys. Chem. A 2015, 119, 5865−5882. (58) Peters, B. Using the Histogram Test to Quantify Reaction Coordinate Error. J. Chem. Phys. 2006, 125, No. 241101.

5138

DOI: 10.1021/acs.jpcb.9b01989 J. Phys. Chem. B 2019, 123, 5131−5138