Multiscale Treatment for the Molecular Mechanism of a DielsâAlder

Article pubs.acs.org/JCTC

Multiscale Treatment for the Molecular Mechanism of a Diels−Alder Reaction in Solution: A QM/MM-MD Study Jorge Soto-Delgado,*,† Ricardo A. Tapia,‡ and Juan Torras*,§ †

Departamento de Ciencias Químicas, Facultad de Ciencias Exactas, Universidad Andrés Bello, Quillota 980, Viña del Mar, Chile Departamento de Química Orgánica, Pontificia Universidad Católica de Chile, Vicuña Mackenna, 4860 Santiago, Chile § Departament d’Enginyeria Química, Escola d’Enginyeria de Barcelona Est (EEBE), Universitat Politècnica de Catalunya, C. Eduard Maristany 10-14, 08019 Barcelona, Spain ‡

S Supporting Information *

ABSTRACT: Thermodynamics and the solvent role in the acceleration of the Diels−Alder reaction between cyclopentadiene (CPD) and methyl vinyl ketone (MVK) have been revisited. In this work we use an ab initio hybrid QM/ MM-MD scheme combined with multiple steered molecular dynamics to extract the free energy pofile in water and methanol using the bidirectional Minh−Adib estimator. We obtain 18.7 kcal mol−1 and 20.8 kcal mol−1 free energy barrier for the reaction in water and methanol, respectively. This methodology reproduces experimental values with an absolute error of about 0.8 kcal mol−1. The experimental difference between the activation free-energy barriers of water and methanol is also reproduced with an absolute error of about 0.1 kcal mol−1. We explore the charge transfer evolution along reaction coordinates to characterize the electronic behavior for this reaction. It is shown that the solvent molecules around the reaction system produce a global polarization along the reaction coordinate which is consistent with the solvent polarity. The results highlight the role of hydrogen bonding formed in the transition state to stabilize the system charge reorganization in the reaction process.

■

INTRODUCTION The Diels−Alder (DA) reaction is the most convenient strategy used in organic synthesis to obtain unsaturated six-membered rings.1 The mechanism of this reaction has been the subject of strong research, and it is suggested that the reaction is concerted; however, stepwise mechanisms are also proposed.2 Consequently, the mechanism as well as solvent effects has been studied using several computational methodologies to explain experimental results. Experimental data suggest that responsible factors on the acceleration of the DA process involve enhanced hydrogen bonding between the solvent and the transition state (TS) structure relative to the initial state of the reaction.3 Moreover, during the activation process, the enforced hydrophobic interactions decrease on the reactant hydrophobic surface area. In this context, the reaction rates for the DA reaction between cyclopentadiene (CPD) and methyl vinyl ketone (MVK) in different solvents have been reported.4 The acceleration of this reaction in water was first obtained by Breslow, who found that both reactivity and endo/exo selectivity are enhanced by water.4 Thus, this reaction resulted in an appreciable decrease in the reaction rate using methanol as solvent (see Scheme 1). Previous theoretical investigations showed the separation and quantification of the relative terms that contribute to the enhanced reaction rates using an explicit treatment of solvent © 2016 American Chemical Society

Scheme 1. Detail of Diels−Alder Reaction and Experimental Data Reported

under the QM/MM framework.5 Jorgensen and co-workers studied the environmental impact of solvent on the reaction by using a combined PDDG/PM3/MM/Monte Carlo method, but theoretical results overestimate the experimental activation barrier height by ∼10 kcal/mol.5c However, the results reported using this methodology revealed a stronger hydrogen bond (HB) for the TS in water when compared with those observed for the reactant and product, mainly due to diene to dienophile charge transfer and a coexistent polarization of the carbonyl groups at the TS.5c Recently, Houk and co-workers6 describe Received: August 3, 2016 Published: September 17, 2016 4735

DOI: 10.1021/acs.jctc.6b00772 J. Chem. Theory Comput. 2016, 12, 4735−4742

Article

Journal of Chemical Theory and Computation

were used. Finally, the system was equilibrated up to a constant density for an additional 500 ps at constant pressure (1 atm) with a pressure relaxation time of 2 ps and a cutoff of 10 Å. All of these steps used the SHAKE algorithm to restrain bonds containing hydrogen14 and an additional restraint for the solute (CPD + MVK) to keep constant an initial distance of 4 Å between CPD and MVK fragments. The reaction coordinate variable (ξ) was built by considering a lineal combination of d1 and d2 distances between carbon atoms of CPD and MVK fragments involved in the reaction process (Figure 1). More specifically, ξ = c1d1 + c2d2, where d1 and d 2 represent the C1−C2 and C3−C4 distances, respectively.

the reaction process using the solvent-perturbed transition state (SPTS) sampling scheme for simulating chemical reaction dynamics in the condensed phase. This method indicated that the HB formation between the solvent and the TS in the reaction between CPD and MVK in water occurs by reorientation of a solvent molecule during the reaction process. In fact, the TS is stabilized by means of a preferential HB, being the largest contribution to the reaction rate acceleration in water. On the other hand, the authors suggest that water polarization surrounding the TS is essential to create a correct representation of dynamical formation of HB in the TS by means of water reorientation.6 Recently, a distinctive multiscale approach has been applied to the reaction kinetics study in order to obtain more accurate thermodynamics and kinetic data of the process using higher ab initio treatment of the system. The hybrid quantum mechanical/molecular mechanical molecular dynamics (QM/ MM-MD) approach has been previously used to study a reaction mechanism in proteins and related biological compounds.7,8 In this multiscale approach the QM forces on atoms are computed in a small subregion of the system using ab initio methodology and the remaining molecules in the MM region are calculated using classical potentials. Finally, the two systems and the QM/MM coupling forces are integrated along simulation time to obtain a dynamics trajectory. Within this context, combining QM/MM-MD and Multiple Steered Molecular Dynamics (MSMD) was used to extract the free energy reaction profile of the process by means of a Jarzynski relationship.9 In this work we will study the DA reaction between cyclopentadiene (CPD) and methyl vinyl ketone (MVK) in water and methanol (Scheme 1) in order to provide a deep insight into the thermodynamics and a clear dynamical picture of the HB role along the reaction coordinate (RC). We have organized the paper in the following way: First, a detailed study using QM/MM-MD methodology, with a broad sampling and incorporating explicit solvent effects, will be conducted to obtain accurate ab initio free energy reaction profiles. Second, a dynamical study on the TS bond formation will be performed by means of a time-evolution analysis of the charge fragments. Finally, the role that the HB plays to stabilize the transition state (TS) will be discussed.

Figure 1. Detail of carbon atom distances d1 and d2 involved in the DA reaction coordinate between methyl vinyl ketone (MVK) and cyclopentadiene (CPD).

Generation of an Equilibrated Initial Ensemble. The Jarzynski equality assumes that all pulling processes start from structures belonging to the same initial equilibrated ensemble. Initial structures were obtained by restraining d1 and d2 distances close to the equilibrium (req) by applying a harmonic potential of 1500 kcal mol−1 Å−2 both centered at 4.0 Å (d1 and d2, respectively). Thus, starting from the previous equilibrated systems in methanol and water, a final production trajectory was obtained to derive an initial set of independent snapshots as a precursor of MSMD along ξ. Classical MD was conducted in an NVT ensemble with periodic conditions at 298 K during 10 ns. Twenty different snapshots were extracted from the whole classical production trajectory as initial structures. To generate free energy reaction profiles (FEP), QM/MM-MD calculations were conducted using an Amber/PUPIL/Gaussian interface15 with the General Amber Force Field (GAFF)10 describing the classical (MM) region and both MPWB1K16 and B3LYP17 functionals together with the 6-31G(d) basis set, as implemented in Gaussian 03,13 to describe the quantum region (QM). After the classical equilibration, the solvated system was divided in two subsystems treated with different levels of calculations. Indeed, the MVK and CPD molecules were treated as quantum atoms (QM region), and the remainder solvent molecules were treated classically (MM region). Initially, an equilibration of 1000 steps within the QM/MMMD approach was performed so that the starting structures were consistent with this level of theory. Multiple Steered Molecular Dynamics (MSMD), Forward Reaction. The FEP is obtained by statistical treatment of the work function along ξ using several trajectories via the MSMD method, which relates the system’s nonequilibrium dynamics to its equilibrium properties.18 MSMD considers an external time-dependent perturbation applied to the system [λ= λ(t)]. The free energy changes ΔG(λ) and the external work performed on the system W(λ) are defined as it evolves from

■

COMPUTATIONAL DETAILS Parameterization and System Preparation. All parameters were taken, with the exception of the atomic partial charges, from the General Amber Force Field (GAFF)10 implemented in the standard Amber 12 distribution.11 The CPD and MVK atomic charges were adjusted through the RESP (restrained electro-static potential) strategy,12 in which the electrostatic potential surface is computed by quantum mechanics methods and fitted using atom-centered point charges. The geometry of each system was optimized at HF/ 6-31G(d) level of theory using the Gaussian 03 suite of programs.13 The system was comprised of the reaction system CPD + MVK and a total of 1760 water or 888 methanol explicit molecules under periodic boundary conditions. Initially the system was prepared using classical MD. Thus, the system was first subjected to 15 000 minimization steps, using the MM Hamiltonian to remove any clashes, and then heated from 0 to 298 K during 60 ps in an NVT ensemble. A time step of 1.0 fs and a Langevin thermostat with a 1.0 ps−1 collision frequency 4736


Article


Figure 2. Free energy reaction profiles of Diels−Alder reaction at MPW1BK/6-31G(d) level in (a) water and (b) methanol. Unidirectional Hummer−Szabo estimator is applied to forward (FW, black line) and reverse (RV, red dotted line) pulling. The bidirectional Minh−Adib estimator (blue line) are also shown.

an initial to a final state (λ0 → λ). Thus, the Jarzynski relationship9 and the Crooks fluctuation theorem19 allow equilibrium free-energy differences to be related to the nonequilibrium work distribution, which, indeed, allows recovering the arbitrary equilibrium ensemble averages from measurements of driven nonequilibrium processes. Based on those theorems there are two different Potential of Mean Force (PMF) estimators that will be applied in this work: the unidirectional Hummer−Szabo estimator20 and the bidirectional PMF estimator of Minh and Adib21 that optimally includes trajectories from the forward and reverse perturbation along ξ. e−βG0(λ) = ⎡ n δ(λ − λ )e−βW0t t ∑t ⎢ F ⎢⎣ nF + nR e−β(W −ΔF)

+

nR δ(λ − λt )e−βW0

nF + nR e F β[V (λ ; t ) −ΔFt ]

t

β(W +ΔF )

thermalized using both DFT functionals in the QM region for 1000 steps (1 ps) and then pulled in the direction of the reverse reaction, i.e., to broken d1 and d2 atomic bonds, following the same procedure as described above. Radial Distribution Functions. ξ values of about ξ = 2.0, 1.2, and 0.8 were chosen as representatives of the reactants, transition state, and products, respectively. Specifically, a global set of snapshots around those ξ values from each of the 20 trajectories were combined to yield a total of 11 ps MD trajectory for the calculation of the Radial Distribution Functions (RDFs), atomic Mulliken charges, and averaged geometries. Non-Covalent Interactions. Characterization of the weak Non-Covalent Interactions (NCI) on the temporally averaged complex structures was performed with the NCIPlot program.24 The NCI surface enables the study of domains of the electronic density associated with weak interactions, being able to distinguish the strength and the attractive or repulsive nature of such interactions.

⎤ ⎥e β ΔFt ⎥ R⎦

∑t e

■

(1)

RESULT AND DISCUSSION The present study has been divided into three parts: (i) First, a complete characterization of the energy profile for the molecular mechanism of the CPD and MVK in methanol and water using B3LYP and the MPW1BK functional are performed; (ii) second, the charge transfer process is evaluated as the sum of the Mulliken charges for each fragment (CPD and MVK) and variation of the local charges along reaction coordinate; and (iii) finally, an analysis of the geometries and solvent effect through radial distribution function and NCI isosurface for TSs are performed in order to explain the increase on the reaction rates experimentally observed. Energetic Aspects. The free energy reaction paths for both solvents are shown in Figure 2 (MPWB1K) and in Figure S1 (B3LYP, in the Supporting Information). Both PMF estimators, Hummer−Szabo and Minh−Adib, were used to reconstruct the unbiased PMF along the ξ for each one of the unidirectional profiles using the first estimator and optimally combine the forward and backward profiles to accelerate the convergence of the free energies on the last estimator. The unidirectional Hummer−Szabo estimator presents a clear biased behavior. The further away from the starting equilibrium state, the more biased the unidirectional PMF, showing a clear overestimation of PMF due to an error accumulation through the trajectory. On the other hand, the free-energy differences

where ΔFt = Ft − F0 is the free energy difference between the equilibrium states defined at time t and 0, respectively. V(λ;t) is the time-dependent potential acting along the reaction coordinate. The Minh−Adib PMF estimator has been shown able to significantly reduce bias in unperturbed free-energy profiles, being previously used in several works such as stretching experiments and binding/unbinding processes.22 However, the use of pulling simulations to accurately reconstruct unperturbed FEP associated with a finite number of nonequilibrium trajectories, under relatively fast pulling speeds, is still a matter of study.23 Starting with the structures properly equilibrated close to the equilibrium distance, d1 and d2 bonds were formed by running independent QM/MM-MD trajectories starting from each of the 20 initial structures obtained as described above, treating the CPD and MVK at both level of theory. In each simulation, a harmonic potential with a large force constant (1500 mol−1 Å−2) was applied, and the center of the potential is moved from 4.0 to1.5 Å at ∼0.46 Å ps−1 (5460 MD steps of 1.0 fs). Reverse Reaction. One of the 20 final structures, after the forward reaction process obtained by pulling from products to reactants, was chosen. The system is then restrained at the final distance 1.45 Å (d1 and d2 distances) and re-equilibrated for 10 ns in an NVT ensemble using classical MD to extract 20 new independent structures. These new initial structures were again 4737


Article


Figure 3. Free energy reaction profiles of Diels−Alder reaction derived using Minh−Adib estimator at (a) MPW1BK/6-31G(d) and (b) B3LYP/631G(d) levels of theory in water (blue) and methanol (red).

barrier obtained at the MPWB1K level is the closest to experimental values in both solvents (see Table 1), also having reproduced its experimental solvent differences, ΔΔG‡Exp = −2.0. More specifically, when the MPWB1K value is compared with the experimental value the differences are 0.9 kcal mol−1 and 0.8 kcal mol−1, whereas at the B3LYP level the free energy barrier differences are 2.1 kcal mol−1 and −0.8 kcal mol−1, for water and methanol, respectively. The good results obtained using MPWB1K might be explained due to a better representation of stacking interactions and hydrogen bonding that is presented in this functional.16,25 On the other hand, FEP obtained using the MPWB1K DFT functional shows an exothermic reaction value of ΔGreac = −21.6 kcal mol−1 and −22.3 kcal mol−1 for methanol and water, respectively, within the expected range of values. B3LYP leads to worse results. From data in Table 1, it is observed that B3LYP emphasizes the difference in free activation energies between both solvents whereas higher values of thermodynamic reaction energies were obtained when compared with those derived from MPW1K. Similarly, it was reported that theoretical calculations of similar cycloadditions, using B3LYP methodology, lead to an overestimation on the reaction exothermicities when comparing with experimental data.26 Charge Transfer Process. The charge transfers along ξ have been represented as the sum of atomic Mulliken charges for each fragment, as a function of reaction coordinate (ξ), and averaged in all simulations (forward and reverse pulling). Figure 4 (MPWB1K results) and Figure S4 (B3LYP, Supporting

between reactant, transition, and product complexes using the Minh−Adib estimator do not present the previous observed biasing on the PMF profile because of the use of both forward and backward reaction profiles. The pulling started from reactants at 4 Å (ξ = 4) and stopped at 1.45 Å (ξ = 0.75) for the forward reaction. Figure 3 shows the FEP of the DA reaction in both solvents and using different DFT functionals. All the work functions for the forward and reverse reaction are also shown in Figure S2 and Figure S3 (Supporting Information). It is observed that at both of the ab initio levels of theory FEP in water presents a lower energy barrier than the FEP in methanol, as expected from experimental data. However, a close inspection of energy values listed in Table 1 allows us to better quantify the Table 1. Free Energy Differences (ΔG, kcal mol−1) Obtained from Multiple Steered Molecular Dynamics (MSMD) and Derived Using the Bidirectional Minh−Adib Estimator method

solvent

ΔG‡Calc

ΔG‡Exp

ΔGreac

MPWB1K

water methanol water methanol

18.7 20.8 17.5 22.4

19.6 21.6 19.6 21.6

−27.6 −26.7 −13.0 −11.8

B3LYP

difference in the energy barrier. The activation free energy barrier differences between methanol and water are ΔΔG‡Calc = −2.1 and −4.9 kcal mol−1 at MPWB1K and B3LYP levels of theory, respectively. Interestingly, the activation free-energy

Figure 4. Charge evolution of fragments along reaction coordinates for DA reaction (a) in water and (b) in methanol derived at MPWB1K/631G(d) level of theory. CPD (black line) and MVK (red line) fragments are shown. 4738


Article


Figure 5. Charge evolution of atoms involved in the process along RC(ξ) for the DA reaction (a) in water and (c) in methanol derived at the MPWB1K/6-31G(d) level of theory. Transition state region is highlighted in a detailed graph for (b) water and (d) methanol.

Information) show the charge evolution along ξ for each one of the fragments. At the CPD−MVK equilibrium distance (ξ = 2) the charge for each moiety is about the same (∼−0.3 e) using the MPWB1K functional. As the C1−C2 and C3−C4 bond distances decrease, there is initially a slight polarization in the reactant direction. Thus, in the TS region the polarization produced a change in the polarization charge transfer between fragments with values of ∼0.12 e and ∼0.18 e for methanol and water, respectively. After the TSs most of the charges in the products at each fragment are stabilized and the bond can be considered formed. This behavior explains the large acceleration observed in these polar DA reactions.27 This also shows that the polarization of reagents is slightly higher in water, which is consistent with the higher polarity of the reaction medium. From the electronic point of view of the different atoms involved in the bond formation process, in Figure 5 we have presented the averaged charges on the carbon atoms of CPD and MVK fragments along ξ in solution at the MPWB1K/631G(d) level (Figure S5 in Supporting Information shows B3LYP results). It can be observed that the charge of all atoms is almost invariant along ξ from the beginning of the reaction. We would like to note, however, about four different points on the transition state region, denoted as P1, P2, P3, and P4 (Figure 5). In the transition state region denoted as P1 that the reaction process begins with an increment on the negative charge of the C1 atom is observed. Consequently, the charge of the C2 atom decreases, leading an initial polarization along C1−C2. After P1, the C1 atom increases its local charge along the first stage of the reaction as a precursor of an initial strong interaction that could be attributed to the first C1−C2 σ bond. Notice that, together with the increment in the negative charge of the C1 atom, a decrease in the local charge for C4 (point P2) takes

place. So, in the transition state region denoted as P1 and P2 the electronic reorganization of the CPD fragment takes place. Consequently, this process produces a large polarization in the C2−C3 bond characterized by the increment in the negative charge of the C3 atom (MVK fragment, point P3). In fact, in the course of the reaction path a tendency to increment the negative charge of the C3 atom is observed; meanwhile, the C4 atom decreases its negative charge. This charge reorganization could be the consequence of a back-donation along the formation of the second C3−C4 σ bond. In addition, this process produces a reduction in the negative charge of C5 (carbonyl carbon) corresponding to point P4 in Figure 5. However, we want to note the observed charge variation in the oxygen atom of the carbonyl group of the MVK fragment, depending on the solvent media. There is a more pronounced increase in the negative charge of the O atom (Figure 5a) in water solvent when compared with the charge evolution of oxygen atom in methanol (Figure 5b). These results are consistent with an HB presence formed during the TS, as reflected in the literature, within several cycloaddition processes.28 HB formed during the TS between the carbonyl oxygen atom and the water hydrogen atoms helps to polarize the CO bond and stabilize the system charge reorganization. Thus, the microsolvation effect using explicit water explains the decrease in the activation barriers due to an electronic polarization of the carbonyl group at the TS. Similarly, the effect of protic solvents on the rate of DA reactions has been explained through HB stabilization of the endo transition state and reduction of the hydrophobic surface area.29 Therefore, near the TS region most of the DA reactions that present a higher asynchronous bond formation process can be considered as a one-step two-stage mechanism, which is characterized by zwitterionic species that are stabilized in polar solvents.27 The global process is detailed in the Scheme 2. 4739


Article

Journal of Chemical Theory and Computation Scheme 2. Bond Formation Evolution in the Transition State Region

Geometric Aspects and Solvent Interactions. The extent of the asynchronicity associated with the bond-formation can be measured by means of the difference between the lengths of the two new bonds that are being formed in the reaction, i.e., Δd = d2(C3−C4) − d1(C1−C2) at the TSs. Considering the averaged bond length obtained in the TSs (Table S1 in Supporting Information), the Δd values obtained on the forward reaction path at the MPWB1K level are 0.45 and 0.41 in methanol and water, respectively. Moreover, similar values of 0.52 and 0.56 are obtained at the B3LYP level in methanol and water, respectively. These results indicate TSs with a highly asynchronous bond-formation process, consistent with that observed in the charge transfer analysis. Figure 6

shows the different averaged bond distance values for the new bonds formed at the TS in methanol and water. From the figure, one can easily observe the asynchronicity on the new bonds formed and HBs that are able to stabilize the TS between the solvent and the carbonyl oxygen of the MVK fragment. Figure 7 shows the radial distribution functions (RDF) of methanol and water molecules around the carbonyl oxygen of the MVK fragment for the reactants, transition state, and product. The small change in the RDFs from the reactants to transition states is mainly due to the reduction of charge in those atoms. Since the transition state structure is close to the products, one would expect the products and TS RDFs to be mostly similar, and this is clearly what is observed in both cases (Table S2, Supporting Information). There is a much stronger HB in water when it is compared with methanol, with maximum distance probabilities of 1.78 and 2.55 Å, respectively, that would help to stabilize the TS in water. Moreover, it is observed from the running number (Table S2) that the number of hydrogen atoms interacting with the CO oxygen in the TSs are 2.23 in water and 2.27 in methanol, within the first solvation cell. However, meanwhile methanol keeps its coordination number along ξ stable, and the DA reaction in water shows different hydrogen coordination numbers of 2.08, 2.23, and 2.26 for reactants, TS, and products, respectively. Thus, there is an increase in the HB occurrence from the TS. Moreover, the HB interactions in water are longer lived than in methanol. In fact, the analysis of HB occupation (considering geometric criteria of distance and angle) through the TS region shows a prominent difference in percentage when both solvents are compared; see Table 2.

Figure 6. Representative average structures for TSs of DA reaction (a) in water and (b) in methanol. Averaged bond distances and its standard deviation are also shown.

Figure 7. Radial distribution function (RDF) between the hydrogen atom and the carbonyl oxygen atom of the MVK fragment in (a) water solvent and (b) methanol solvent. Reactants, TS, and product RDFs are also shown. 4740


Article


better results with less biased profiles using fast velocity pulling to obtain ab initio reaction profiles. The MPWB1K DFT functional gives a better representation of the kinetics and thermodynamics reaction parameters, with free-energy barriers ΔG‡ = 18.7 and 20.8 kcal mol−1 for water and methanol, respectively, very close to experimental values with an absolute error of 0.9 and 0.8 kcal mol−1 in water and methanol, respectively. The acceleration of the DA reaction in water is also observed with an activation free-energy barrier difference, ΔΔG‡ = 2.1 kcal mol−1 between methanol and water, with an error of 0.1 kcal mol−1. The electronic behavior of the DA reaction can be explained through a global polarization along reaction coordinates produced by the solvent molecules being consistent with the polarity of the reaction medium. The charge evolution of the atoms involved in the reaction process indicates that the DA reaction proceeds in a one-step two-stage mechanism through asynchronous TS. The stabilization of zwitterionic species formed in the TS is favored in the polar solvent, whereas asynchronicity of distances in the TS bond formation has similar values in both solvents. Consequently, the HB formed during the TS between the carbonyl oxygen atom and hydrogen atoms in solvent is responsible to polarize the C O bond and stabilizes the system charge reorganization in the reaction process.

Table 2. Averaged Occupation (%) of Hydrogen Bonding between Hydrogen Atom and Carbonyl Oxygen Atom of MKV Fragment through All the Trajectories in the Transition State Region method

solvent

HBa

HBb

HBc

MPWB1K

water methanol water methanol

92.6 94.6 52.1 54.8

60.7 65.9 -

17.6 7.9 -

B3LYP

Weak noncovalent interactions (NCI) between the averaged structures for TSs in both solvents were analyzed by examining the reduced electron density gradient with the NCIPlot program (Figure 8).24 This methodology provides a useful

■

Figure 8. NCI plots for Solute−Solvent interaction (a) in water and (b) methanol for average structures.

ASSOCIATED CONTENT

S Supporting Information *

description of molecular interactions, which are frequently represented by an arbitrary color code: blue, green, and red, which are used for highly attractive weak interactions (such as hydrogen bonds), extremely weak interactions (such as van der Waals), and repulsive interactions (such as steric clashes), respectively. NCI mapping figures allows an easy identification of the regions with strong and weak electron pairing. Figure 8 displays the reduced density gradient isosurface of the TS region in methanol and water, respectively. Both TSs are dominated by van der Waals interactions between the solvent and the TSs−complex. Note that the solvent−complex interactions are also dominated by nonspecific van der Waals interactions in the two studied environments. However, the complex in water presents three localized blue regions at the isosurfaces in comparison with two localized blue regions in methanol. These results suggest that the formation of specific HBs is more effective in water with respect to methanol. This is also consistent with the observed number of coordinated hydrogen atoms with the carbonyl oxygen, shown above. The role the weak interactions plays in the stability of the complex is higher in water than in methanol, which also indicates that the factors responsible for the acceleration of the DA process involves enhancing hydrogen bonding between the solvent and the TS structure relative to the initial state of the reaction.

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00772. Tables of averaged bond distances and coordination geometries of involved species and figures of the work and PMF reaction profiles of the system reactivity in water and methanol (PDF)

■

AUTHOR INFORMATION

Corresponding Authors

*(J.S.-D.) E-mail: [email protected]. *(J.T.) E-mail: [email protected]. Funding

J.S.-D. would like to thank the support from FONDECYT Grants 3130359 and 11150988. R.A.T. would like to thank Fondo de Innovación para la Competitividad del Ministerio de Economia,́ Fomento y Turismo, Chile (Project RC 130006 CILIS). J.T. would like to thank MINECO-FEDER (MAT2015-69367-R) and DIUE of the Generalitat de Catalunya (XRQTC and CESCA) for funding. Notes

The authors declare no competing financial interest.

■ ■

■

ACKNOWLEDGMENTS J.T. acknowledges PRACE for access to Curie TN and HN, France, at GENCI@CEA and Mare Nostrum, Spain.

CONCLUSIONS The molecular mechanism of the DA reaction between cyclopentadiene (CPD) and methyl vinyl ketone (MVK) has been studied at the B3LYP/6-31G(d) and MPWB1K/631G(d) levels of theory using a QM/MM-MD combined with MSMD approach in a set of forward and reverse bond pulling along the reaction coordinates. A free-energy reaction profile was derived in water and methanol using unidirectional Hummer−Szabo and bidirectional Minh−Adib Potential of Mean Force estimators. The Minh−Adib estimator presents

REFERENCES

(1) Carruthers, W. Cycloaddition Reactions in Organic Synthesis; Pergamon Press: Oxford, U.K., 1990. (2) (a) Domingo, L. R. RSC Adv. 2014, 4, 32415−32428. (b) Yepes, D.; Murray, J. S.; Pérez, P.; Domingo, L. R.; Politzer, P.; Jaque, P. Phys. Chem. Chem. Phys. 2014, 16, 6726−6734. (c) Loncharich, R. J.; Brown, F. K.; Houk, K. N. J. Org. Chem. 1989, 54, 1129−1134. (d) Birney, D. M.; Houk, K. N. J. Am. Chem. Soc. 1990, 112, 4127−4133. 4741


Article

Journal of Chemical Theory and Computation (e) Sustmann, R.; Sicking, W. J. Am. Chem. Soc. 1996, 118, 12562− 12571. (3) (a) Blokzijl, W.; Blandamer, M. J.; Engberts, J. B. F. N. J. Am. Chem. Soc. 1991, 113, 4241−4246. (b) Otto, S.; Blokzijl, W.; Engberts, J. B. F. N. J. Org. Chem. 1994, 59, 5372−5376. (c) Wijnen, J. W.; Engberts, J. B. F. N. J. Org. Chem. 1997, 62, 2039−2044. (4) (a) Rideout, D. C.; Breslow, R. J. Am. Chem. Soc. 1980, 102, 7816. (b) Breslow, R.; Guo, T. J. Am. Chem. Soc. 1988, 110, 5613. (c) Breslow, R. Acc. Chem. Res. 1991, 24, 159−164. (5) (a) Blake, J. F.; Lim, D.; Jorgensen, W. L. J. Org. Chem. 1994, 59, 803−805. (b) Chandrasekhar, J.; Shariffskul, S.; Jorgensen, W. J. J. Phys. Chem. B 2002, 106, 8078−8085. (c) Acevedo, O.; Jorgensen, W. L. J. Chem. Theory Comput. 2007, 3, 1412−1419. (6) (a) Yang, Z.; Doubleday, C.; Houk, K. N. J. Chem. Theory Comput. 2015, 11, 5606−5612. (b) Liu, F.; Yang, Z.; Mei, Y.; Houk, K. N. J. Phys. Chem. B 2016, 120, 6250−6254. (7) (a) Xiong, H.; Crespo, A.; Marti, M.; Estrin, D.; Roitberg, A. E. Theor. Chem. Acc. 2006, 116, 338−346. (b) Lonsdale, R.; Harvey, J. N.; Mulholland, A. J. Chem. Soc. Rev. 2012, 41, 3025−3038. (c) Ramírez, C. L.; Zeida, A.; Jara, G. E.; Roitberg, A. E.; Martí, M. A. J. Chem. Theory Comput. 2014, 10, 4609−4617. (d) Sherwood, P.; Brooks, B. R.; Sansom, M. S. P. Curr. Opin. Struct. Biol. 2008, 18, 630−640. (8) (a) Torras, J.; Roberts, B. P.; Seabra, G. M.; Trickey, S. B. Adv. Protein Chem. Struct. Biol. 2015, 100, 1−31. (b) Torras, J.; Seabra, G. M.; Roitberg, A. E. J. Chem. Theory Comput. 2009, 5, 37−46. (9) Jarzynski, C. Phys. Rev. Lett. 1997, 78, 2690−2693. (10) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157−1174. (11) Case, D. A.; Darden, T. A.; Cheatham, T. E., III; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M.; Roberts, B.; Hayik, S.; Roitberg, A.; Seabra, G.; Swails, J.; Götz, A. W.; Kolossváry, I.; Wong, K. F.; Paesani, F.; Vanicek, J.; Wolf, R. M.; Liu, J.; Wu, X.; Brozell, S. R.; Steinbrecher, T.; Gohlke, H.; Cai, Q.; Ye, X.; Wang, J.; Hsieh, M.-J.; Cui, G.; Roe, D. R.; Mathews, D. H.; Seetin, M. G.; Salomon-Ferrer, R.; Sagui, C.; Babin, V.; Luchko, T.; Gusarov, S.; Kovalenko, A.; Kollman, P. A. AMBER 12; University of California: San Francisco, 2012. (12) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269−10280. (13) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, Revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (14) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327−341. (15) (a) Torras, J.; Deumens, E.; Trickey, S. B. J. Comput.-Aided Mater. Des. 2006, 13, 201−212. (b) Torras, J.; Seabra, G. M.; Deumens, E.; Trickey, S. B.; Roitberg, A. E. J. Comput. Chem. 2008, 29, 1564−1573. (16) Zhao, Y.; Truhlar, D. G. J. Phys. Chem. A 2004, 108, 6908−6918. (17) (a) Becke, A. D. J. Chem. Phys. 1993, 98, 5648−5652. (b) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785−789.

(18) (a) Isralewitz, B.; Gao, M.; Schulten, K. Curr. Opin. Struct. Biol. 2001, 11, 224−230. (b) Jensen, M. Ø.; Park, S.; Tajkhorshid, E.; Schulten, K. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 6731−6736. (19) Crooks, G. E. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1999, 60, 2721−2726. (20) (a) Hummer, G.; Szabo, A. Proc. Natl. Acad. Sci. U. S. A. 2001, 98, 3658−3661. (b) Hummer, G.; Szabo, A. Acc. Chem. Res. 2005, 38, 504−513. (21) Minh, D. D. L.; Adib, A. B. Phys. Rev. Lett. 2008, 100, 180602− 180604. (22) (a) Do, T. N.; Carloni, P.; Varani, G.; Bussi, G. J. Chem. Theory Comput. 2013, 9, 1720−1730. (b) Frey, E. W.; Li, J.; Wijeratne, S. S.; Kiang, C.-H. J. Phys. Chem. B 2015, 119, 5132−5135. (23) Ngo, V. A.; Kim, I.; Allen, T. W.; Noskov, S. Y. J. Chem. Theory Comput. 2016, 12, 1000−1010. (24) (a) Johnson, E. R.; Keinan, S.; Mori-Sánchez, P.; ContrerasGarcía, J.; Cohen, A. J.; Yang, W. J. Am. Chem. Soc. 2010, 132, 6498− 6506. (b) Contreras-García, J.; Johnson, E. R.; Keinan, S.; Chaudret, R.; Piquemal, J.-P.; Beratan, D. N.; Yang, W. J. Chem. Theory Comput. 2011, 7, 625−632. (25) Dkhissi, A.; Blossey, R. Chem. Phys. Lett. 2007, 439, 35−39. (26) (a) Guner, V.; Khuong, K. S.; Leach, A. G.; Lee, P. S.; Bartberger, M. D.; Houk, K. N. J. Phys. Chem. A 2003, 107, 11445− 11459. (b) Pieniazek, S. N.; Clemente, F. R.; Houk, K. N. Angew. Chem., Int. Ed. 2008, 47, 7746−7749. (27) Domingo, L. R.; Sáez, J. A. Org. Biomol. Chem. 2009, 7, 3576− 3583. (28) (a) Soto-Delgado, J.; Aizman, A.; Contreras, R.; Domingo, L. R. Molecules 2012, 17, 13687−13703. (b) Frapper, G.; Bachmann, C.; Gu, Y.; Coval De Sousa, R.; Jérôme, F. Phys. Chem. Chem. Phys. 2011, 13, 628−636. (c) Gordillo, R.; Dudding, T.; Anderson, C. D.; Houk, K. N. Org. Lett. 2007, 9, 501−503. (d) Polo, V.; Domingo, L. R.; Andrés, J. J. Phys. Chem. A 2005, 109, 10438−10444. (29) (a) Furlani, T. R.; Gao, J. J. Org. Chem. 1996, 61, 5492−5497. (b) Engberts, J. B.F.N. Pure Appl. Chem. 1995, 67, 823−828. (c) Kong, S.; Evanseck, J. J. Am. Chem. Soc. 2000, 122, 10418−10427.

4742


Multiscale Treatment for the Molecular Mechanism of a DielsâAlder

Recommend Documents

Multiscale Treatment for the Molecular Mechanism of a DielsâAlder