Combined Molecular Dynamics and Coordinate Driving Method for

Subscriber access provided by University of Sunderland

Reaction Mechanisms

Combined Molecular Dynamics and Coordinate Driving Method for Automatic Reaction Pathway Search of Reactions in Solution Manyi Yang, Lijiang Yang, Guoqiang Wang, Yanzi Zhou, Daiqian Xie, and Shuhua Li J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.8b00799 • Publication Date (Web): 23 Oct 2018 Downloaded from http://pubs.acs.org on October 25, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Combined Molecular Dynamics and Coordinate Driving Method for Automatic Reaction Pathway Search of Reactions in Solution Manyi Yang, † Lijiang Yang,‡ Guoqiang Wang,† Yanzi Zhou,*,† Daiqian Xie† and Shuhua Li*,† †School

of Chemistry and Chemical Engineering, Key Laboratory of Mesoscopic

Chemistry of Ministry of Education, Institute of Theoretical and Computational Chemistry, Nanjing University, Nanjing, 210023, People’s Republic of China ‡

Beijing National Laboratory for Molecular Sciences, College of Chemistry and

Molecular Engineering, Peking University, Beijing 100871, People’s Republic of China

KEYWORDS: MD/CD method, automatic, reaction pathway search, reactions in solution ABSTRACT: The combined molecular dynamics and coordinate driving (MD/CD) method is further refined in this work to automatically search reaction pathways for chemical reactions in solution. In this refined MD/CD method, the selective integrated tempering sampling based MD (SITS-MD) simulations are performed to efficiently sample conformers of the reactants in a realistic solution environment. Then, dozens of the reactant/solvent clusters from SITS-MD simulations were used in the subsequent CD step

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

with the quantum mechanics/molecular mechanics (QM/MM) method to generate reaction pathways. The present MD/CD method is able to search reaction pathways, in which solvent molecules may directly participate. This approach is applied to investigate two reactions without any prior knowledge of the reaction mechanism: Cope elimination of amine oxide in aqueous and dimethylsulfoxide solutions, and dehydration of methanediol in aqueous solution. For both reactions, our calculations can locate a large number of lowenergy reaction pathways. For the dehydration reaction in aqueous solution, free energy barriers for several reaction modes located by the MD/CD method have been obtained from potential of mean force calculations. Our results show that the most likely reaction mode is the dehydration of methanediol catalyzed by two water molecules. These two illustrative applications demonstrate that the refined MD/CD method is a promising tool in predicting low-energy reaction pathways for reactions in solution.

INTRODUCTION

Theoretical understanding, prediction and design of chemical reactions are very important applications of theoretical chemistry.1,2 These applications can be facilitated by efficient methods of automatically searching reaction pathways. Finding the stationary points including both local minima and transition state (TS) structures along the reaction pathways on the potential energy surface (PES) is essential for understanding the molecular mechanism of chemical reactions. Thanks to the optimization techniques developed for 2


Page 2 of 36

Page 3 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


locating such stationary points, several automatic methods for searching reaction pathways systematically without any initial guess have been proposed. These methods include global reaction root mapping (GRRM) strategy,3,4 with the anharmonic downward distortion following (ADDF) method5-7 and artificial force induced reaction (AFIR) method8-11 proposed by Morokuma and co-workers, the single-ended growing string method (GSM)12,13 proposed by Zimmerman’s group, the stochastic surface walking based reaction sampling (SSW-RS) method14-16 proposed by Liu’s group, the transition state search using chemical dynamics simulations (TSSCDS) method17,18 proposed by Martínez-Nuñez, and so on.19-25 These methods have been used to investigate a variety of chemical reactions. Very recently, our group proposed the combined molecular dynamics and coordinate driving (MD/CD) method,26 in which the cost-effective MD simulations at the molecular mechanics (MM) or semi-empirical quantum mechanics (QM) level are performed to locate conformational isomers for minimum structures and a refined coordinate driving (CD) technique is used to build reaction pathways. In the CD technique originally proposed by Koča and co-workers,24 a series of constrained minimizations by fixing given internal coordinates (but relaxing all the other coordinates) are carried out to drive the reaction in the direction of the desired product. In our refined CD step, numerous torsion or anglebased reaction coordinates are excluded because the corresponding conformational changes are described by classical MD simulations. As a result, the MD/CD method is computationally much more efficient than other automatic methods, and thus is applicable 3



to relatively large chemical reactions. However, most chemical reactions occur in solution. The interaction between solute and solvent has a significant impact on the mechanism of chemical reactions, especially for those reactions in protic or polar solutions. Polarizable continuum model (PCM)27,28 is one of the most widely used methods to deal with the solvent effect on chemical reactions, in which the solvent is modeled as a continuous dielectric of infinite extent that surrounds a cavity containing the solute molecule. This implicit model can provide satisfactory descriptions on the solvent effect for systems with relative weak reactant-solvent interactions at a relatively low computational cost. Unfortunately, this method is less accurate in describing the effect of highly specific reactant-solvent interactions, such as electrostatic or hydrogen-bond interaction. Moreover, sometimes the solvent molecules may directly participate in the reaction29-32. For such reactions, some of the neighboring solvent molecules around the reactants should be explicitly included in QM calculations. In general, the full QM treatment of the reactants and its neighboring solvent molecules is computationally very demanding. A more realistic model for simplifying the calculation is to adopt a microsolvation model, which consists of the reactants and those neighboring solvent molecules in a few solvation shells. It is well known that the hybrid quantum mechanics/molecules mechanics (QM/MM) method is an appropriate theoretical tool to treat a microsolvation model. With this hybrid approach, the reaction center (the reactants and a few spatially neighboring solvent molecules) is calculated at the QM level, while the 4


Page 4 of 36

Page 5 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


rest of the solvent molecules in a microsolvation model is calculated at the MM level. Recently, by combining with the n-layered integrated molecule orbital and molecular mechanics method (ONIOM, an alternative QM/MM method), the GRRM method was applied to explore reaction pathways of an aldol reaction in an aqueous solution33,34. In their microsolvation model, a single water is manually assigned as a catalyst for the Aldol reaction (between H2CO and H2C=C(H)OH) and 299 H2O molecules was added to mimic the local environment of this reaction. However, for a general reaction, it is difficult for users to give predictions on the number of solvent molecules directly participating in the reaction. Therefore, an efficient automatic method for searching reaction pathways of chemical reactions in solution, without any prior knowledge of the reaction mechanism and the number of solvent molecules directly reacting with the reactants, is desired. In the present study, our focus is to automatically search reaction pathways for chemical reactions in solution based on a number of the reactant/solvent clusters. The MD/CD method is a cost-effective theoretical tool for searching low-energy reaction pathways for relatively large chemical reactions. Our test calculations demonstrated that this method performs well in searching reaction pathways for both intermolecular and intramolecular reactions in the absence of explicit solvent molecules. In this paper, the MD/CD algorithm is further refined to allow reaction pathways of chemical reactions in solution to be automatically explored. To better explore the configurational space of the reactants in solution, we employ the selective integrated tempering sampling (SITS)35-38 5



based MD simulations, instead of standard MD simulations used in our previous study. To speed up calculations on reactant/solvent clusters, the ONIOM method is used to generate reaction pathways in the CD step. The refined MD/CD method is then applied to investigate the Cope elimination of anime oxide in both aqueous and DMSO solutions, and the dehydration of methanediol in aqueous solution. In the second reaction, a few water molecules are shown to take part in the reaction as catalysts. Free energy barriers of this reaction for several reaction modes located by the MD/CD method have been obtained from potential of mean force (PMF) calculations. The refined MD/CD method is demonstrated to be an efficient tool for automatically searching low-energy pathways of chemical reactions in solution without any prior knowledge of the reaction mechanism.

2. METHOD

2.1 MD/CD method In the MD/CD method reported previously, the standard MD method is employed for systematically sampling conformational isomers for the minimum structures and a significantly modified CD technique is used to generate bond-making/breaking structural isomers, respectively. The MD/CD method is implemented with a recursive algorithm and each cycle of this method consists of the following steps: (1) Perform MD simulation with a starting structure to produce a set of representative conformers; 6


Page 6 of 36

Page 7 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(2) Use the CD technique for each conformer generated from step 1 to build possible one-step reaction pathways leading to structural isomers (via bond-making/breaking process). Meanwhile, the intermediates and the corresponding TSs along these pathways are obtained; (3) Check whether the energies of either the intermediates or the TSs along the one-step pathways are higher than searching pathway,

Emin  Ecut ( Emin : lowest-energy point along the current

Ecut : relative energy cutoff). If not, go to step 4;

(4) Check whether the obtained intermediate is a new structural isomer or not (details of geometrical comparison between two structures are presented in the Supporting Information (SI)). If yes, add the obtained new structure to the intermediate list. After completion of this cycle, the new structures in the intermediate list become starting structures of the next cycle. The recursive loop will be terminated until all structures in the intermediate list have traversed the four steps above. Finally, the network of different multistep reaction pathways, including all stationary points (intermediates and TSs), can be automatically built. As demonstrated in our previous work26, a comprehensive conformer list of the minimum structures, which is quite important for automatic reaction pathway search, can be obtained via normal MD simulation at an elevated temperature (such as 600K) for systems in gas phase. However, for reactants in condensed phase, simply raising the simulation temperature for the whole system will make the reactants to stay in gas-phase7



like environment, in which the possible reactant-solvent interaction may be broken. Since the SITS method is an efficient tool for selectively enhancing the sampling of the configuration space of the reactants without destroying the solvent environment, here we adopted the SITS-MD simulations (instead of normal MD simulations) to sample conformers of reactants in a realistic solution environment. The schematic picture for SITSMD simulations is illustrated in Figure 1(a), in which only the atoms shown in ball-andstick type are treated with the enhanced sampling. In order to treat the solvent effects accurately, a large number of neighboring solvent molecules are considered here. In the present MD/CD method, the reactant/solvent cluster which contains the reactants and hundreds of neighboring solvent molecules, as shown in Figure 1(b), is used for CD calculations. Given that chemical reactions generally take place within the reaction center (consisting of the reactants and some nearest-neighboring solvent molecules) and full QM calculations of various reactant/solvent clusters are time consuming, we applied the hybrid QM/MM method in the CD technique to perform calculations for the reactant/solvent clusters. Here the ONIOM method implemented in the Gaussian 0939 program, a widely used QM/MM method for treating explicit solvent effects, was used. The details of the SITS-MD simulation and the ONIOM method are given in the following sections.

8


Page 8 of 36

Page 9 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1. Schematic picture of the MD/CD calculations for searching reaction pathways in solution; (a) the reactant and solvent molecules in a periodic box for SITSMD simulation, atoms treated with the enhanced sampling are shown in ball-and-stick type; (b) the reactant/solvent cluster extracted from the SITS-MD trajectory for ONIOM calculations; Atoms included in the QM region are shown in ball-and-stick type. Atoms included in the MM region are shown in wireframe. The atoms in bluewireframe are freely relaxed, and those in red-wireframe are frozen during ONIOM calculations.

2.2 SITS-MD simulation for sampling conformers of the reactants

The SITS-MD simulation is applied in the MD/CD method to sample conformers of reactants in a realistic solution environment. By introducing a potential surface scaled with temperature, SITS enables us to selectively enhance sampling of the configuration space of the solute molecules (the reactants) but keep the solvent molecules at the target 9



temperature

Page 10 of 36

T0 . The potential surface as a function of temperature is shown in Eq. (1), U eff  Eslv 

1

0

ln  nk e

  k ( Eslu  Eslu slv )

(1)

k

Where U eff is the biased potential energy;

Eslo , Eslv and Eslo slv are the internal energy

of the solute, the internal energy of the solvent, and the interaction energy between the

 0  1/ k BT0 , where k B

solute and solvent, respectively; and

is the Boltzmann constant;

k

nk (obtained through an iterative procedure) are the Boltzmann factor and weighting

factor at temperature

Tk ; The biased force Feff can be calculated via Eq. (2), where F

is the force calculated using the original potential function.

Feff =

U eff r

 n e  



k

k

kU

k

 0  nk e

F

(2)

  kU

k

Finally, the desired distribution function  (r ) at the target temperature calculated via Eq. (3), where

 (r )  eff (r )(r ) e

T0 can be

eff (r ) is calculated from the enhanced sampling simulation.   0 [U ( r ) U eff ( r )]

(3)

After the SITS-MD simulation, a certain number ( N c : a predefined parameter) of reactant/solvent clusters (consisting of the reactants and hundreds of neighboring solvent molecules) are extracted from the SITS-MD production trajectory for subsequent reaction pathway searches. For two-component reactants, an additional force is added to restrain the distance between the center of mass of two reactants within a suitable range in our MD simulations. Before CD calculations, all clusters are fully optimized with the ONIOM 10


Page 11 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


method. The details of ONIOM calculations are given in section 2.3. 2.3 CD technique based on ONIOM calculations for building reaction pathways The procedure of the CD technique combined with the ONIOM method for generating reaction pathways consists of the following three steps. First, for a given reactant/solvent cluster, a target reaction coordinate is selected automatically among the specified target atoms. As discussed in our previous paper,26 only atom-atom pairs (for those target atoms) with the condition 0.95  ij  3.5 are regarded as the target reaction coordinates, since a chemical transformation usually occurs only between spatially close atom-atom pairs. Here

ij  rij / ( Ri  R j ) is defined, with

Ri

as the covalent radius of target atom i , and rij as

the distance between target atoms i and j . Then, a serial of constrained ONIOM optimizations, in which a target reaction coordinate is fixed at specified values while all of the other coordinates are freely optimized, are performed to generated an approximate local potential surface. The CD calculation for this target reaction coordinate will be terminated once the energy point along the approximate reaction pathway is higher than

Ecut  Emin .

With this criterion, lots of unreasonable high-barrier pathways can be excluded. The terminal and maximum energy point along this approximate reaction pathway are taken as the approximate intermediate and TS, respectively. Finally, re-optimization of both the approximate intermediate and TS without any constraints gives the true intermediate and TS structures. After completion of the above three steps, a one-step reaction pathway leading to a new structural isomer via the located TS can be built automatically. Meanwhile, 11



Page 12 of 36

this pathway will be removed and the intermediate will be excluded from the list of reaction intermediates if the energy of the located intermediate (or TS) is higher than

Ecut  Emin .

For the ONIOM calculations, the reactant/solvent cluster is divided into two layers: the QM region and the MM region. In general, the reactants and a few nearest-neighboring solvent molecules are calculated at the QM level and the rest of the solvent molecules are calculated at the MM level. As shown in Figure 1 (b), atoms in ball-and-stick type and wireframe are included in the QM and MM region, respectively. To diminish the effect of distant solvent molecules on energies of

stationary structures along the pathways, we

adopt a further simplification, where the inner shell solvents (about 50 solvent molecules) in the MM region (atoms in blue wireframe) are freely relaxed and the outer shell solvent molecules (atoms in red wireframe) are not allowed to move during geometry optimizations. The use of this simplified model for ONIOM calculations allows the optimization of minimum structures and TSs much less costly. Previous work showed that this simplification has little effect on the structures of located stationary points.40,41 It is worth pointing out that no atoms in solvent molecules are considered as active atoms in the present MD/CD method. With fixing the targeted reaction coordinate (consisting of atoms in the reactants) during constrained optimizations (all of the other coordinates are freely optimized), the neighboring solvent molecules (included in the QM region) can automatically take part in the reaction with the reactants (if necessary). Our calculations demonstrated that this MD/CD method can successfully locate the reaction 12


Page 13 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


pathways, in which one or more solvent molecules are involved in the reaction. Since atoms in the solvent molecules are excluded from the target atom list, the freedom of target reaction coordinates for reactions in solution is comparable to those in corresponding gasphase reactions. Thus, this MD/CD method is a cost-effective approach for studying chemical reactions in solution environment. 2.4 PMF calculations for distinguishing the most likely reaction pathway With the refined MD/CD method, one can locate a number of lower-energy pathways for reactions in solution without any prior knowledge of the reaction mechanisms. Generally, to distinguish the most likely reaction pathway in solution, one can perform potential of mean force (PMF) calculations at finite temperatures for several possible lowenergy pathways to estimate the free energy changes along certain reaction coordinates, in which configurational entropy of the reactants and the solvent can be incorporated. Note that in PMF calculations, equilibrium solvation is assumed during the overall reaction process from reactants to products. This assumption may be not valid for all types of reactions. For those reactions, that equilibrium solvation may not occur between the TS and possible products, the use of free energy profiles may not accurately predict the reaction pathways.42,43 In this work, we have combined MD/CD calculations and PMF calculations to obtain the free energy barriers for chemical reactions in solution. The most likely reaction pathway may be defined as the one with the lowest free energy barrier.

13



2.5 Computation details The MD/CD calculations were performed with the automated design of chemical reaction (ADRC) program developed by us. In the ONIOM calculations, density functional theory (DFT) method with suitable bases sets is used for QM calculations and the GAFF2 force field44 (TIP3P45 model for water molecules) is used for MM calculations. The Gaussian 0939 program is used to do ONIOM calculations and all geometry optimizations (with or without constrains) for all stationary points along different reaction pathways. Vibrational frequencies were calculated for all stationary points, including minimum structures (reactants, intermediates, and products) and TSs to confirm whether the optimized geometry corresponds to a minimum or TS. Intrinsic reaction coordinate (IRC) calculations were performed to verify whether a TS structure connects the correct minimum structures. All ONIOM calculations were performed with a mechanical embedding scheme. All energies reported in this paper were corrected by the zero-point vibrational energy. To diminish the energy difference from the nonreactive center, the energy barrier ΔE used in this paper is defined as the difference of the electronic energy between the reactants optimized from the endpoint of the forward IRC profile and the corresponding TS. The AMBER16 package46 is used for performing MD simulations with restrained electrostatic potential (RESP) charges47,48. In this paper, periodic boundary conditions were applied for the reactants/solvent systems. Only atoms in the reactants are treated with the 14


Page 14 of 36

Page 15 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


enhanced sampling. The details of SITS-MD simulations are described as follows. First, the target system was energy minimized for 6000 steps with the steepest descent method and 4000 steps with the conjugated gradient method. Then the system is heated up to 298 K for 500 ps. The integration time step was taken as 1 fs for all MD simulations. Subsequently, the system was equilibrated to an appropriate volume at NPT (P = 1 atm, T = 298 K) ensemble for 500 ps and followed by 2ns of equilibration at constant volume and temperature of 298 K. After the above equilibration steps, 2ns SITS-MD simulation, where

nk ’s in Eq. (1) were updated every 0.1ps, for determining proper nk series for subsequent production run (use constant

nk ’s) was performed. In the SITS-MD simulation, the

temperature used in Eq. 2 is ranged from 298K to 600K with 200 sampling intervals uniformly. Finally, the production trajectory was accumulated for 12 ns, where the first 2 ns was for equilibration, and the last 10 ns was used to generate configurations for subsequent data analysis. As discussed in our previous work26, the use of a large value of representative reactants/solvent clusters ( N c ) may lead to more comprehensive reaction pathways. In this paper,

N c  80 configurations are extracted from the production

trajectory for subsequent reaction pathway sampling. The PMF calculations were performed here by employing Born−Oppenheimer ab initio QM/MM MD49-53 in conjunction with umbrella sampling.53-55 All our QM/MM simulations have been carried out with modified Q-Chem56 and Amber57 programs. For more computational details, please see the SI. 15



3. RESULTS AND DISCUSSION

In this section, we will apply the MD/CD method to investigate the reaction pathways for two reactions without any prior knowledge of the reaction mechanism: Cope elimination of threo-N, N-Dimethyl-3-phenyl-2-butylamine oxide in aqueous and dimethylsulfoxide (DMSO) solutions, and the dehydration of methanediol in aqueous solution. For Cope elimination reaction, each reactant/solvent cluster for ONIOM calculations consists of the amine oxide and 100 solvent molecules, where the amine oxide and its nearest 5 H2O (or 3 DMSO) molecules are included in the QM region and 95 water molecules are included in the MM region. The inner 45 solvent molecules in the MM region are relaxed during optimization (the remaining 50 solvent molecules are fixed). Molecules included in the QM region are treated at the B97X-D58/6-31G* level. All H2O and DMSO in the MM region are treated with the TIP3P model and the GAFF2 force field, respectively. For the second reaction, each reactant/solvent cluster consists of a methanediol molecule and 120 water molecules. Among them, the methanediol and its neighboring 10 H2O molecules are included in the QM region and others are included in the MM region. The inner 50 water molecules in the MM region are relaxed during optimization. Molecules included in the QM region are calculated at the B3LYP59/6-31G** level. All H2O in the MM region are treated with the TIP3P model. Based on the MD/CD results, PMF calculation are performed for several possible reaction pathways. For the QM/MM MD simulation, the QM region consists of methanediol and a few water molecules, which are 16


Page 16 of 36

Page 17 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


directly involved in this reaction. Molecules included in the QM region are described at the B3LYP/6-31G** level. The remaining solvent water molecules were treated molecular mechanically with the TIP3P water model. 3.1 Cope Elimination Reaction in Aqueous and Dimethylsulfoxide Solutions The Cope elimination of amine oxide is a classical reaction for stereospecific synthesis of olefins.60-63 This reaction is ideal for evaluating the efficiency of the refined MD/CD method because it is extraordinarily sensitive to solvent effects (the protic solvent can provide up to million-fold rate decelerations). Unfortunately, for this reaction, the PCM model seriously underestimates the solvent effects in going from the aprotic solvent to the protic solvent.63 Here, the MD/CD method is applied to study this reaction in both aqueous and dimethylsulfoxide (DMSO) solutions with

Ecut  45 kcal/mol. Due to the high

computational cost, only atoms of C12, H13, C14, H15, N24 and O33 are denoted as target atoms (as shown in structure A1 in Figure 3). To illustrate the sampling efficiency of the SITS-MD algorithm, the distribution of conformational isomers of amine oxide in aqueous solution sampled with the SITS-MD method and the normal MD simulations are compared. Since the diversity of conformational change of amine oxide can be described by different values of H13-O33 distance, we have shown the distribution of different conformers as a function of time in Figure 2 (b), and four typical conformers (S1-S4) of the amine oxide in Figure 2 (a). By comparing the distribution of H13-O33 distances from normal MD (298K and 600K) and 17



SITS-MD (298K) trajectories, we found that the SITS-MD simulation can cover a wide conformational space, which is comparable with the normal MD simulation under 600K. Such comparison indicates that a comprehensive conformer list of the reactants can be sampled using the SITS-MD method. In addition, the radial distribution functions (RDFs) of hydrogen bonds (O-HW) between the oxygen of amine oxide and hydrogens of water from SITS-MD and normal MD simulations are compared in Figure 2(c). The results show that the local water environment around the amine oxide in the SITS-MD simulation is quite similar to that in the normal MD simulation under 298K, which indicates that the solvent environment from the SITS-MD simulation is largely retained at the target temperature (298K). Thus, the SITS-MD method can sample a comprehensive conformer list of the minimum structures for amine oxide in aqueous solution.

18


Page 18 of 36

Page 19 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2. Comparison of the sampling efficiency between normal MD simulations (at 298K and 600K), and SITS-MD simulations (at 298K). (a) Typical conformers of amine oxide with different H13-O33 distances (in Å). (b) Distribution of H13-O33 distances as a function of time; (c) Computed RDFs between the oxygen of amine oxide and hydrogens of water (O-HW).

For this Cope elimination reaction in the aqueous solution, our calculations have produced 36 pathways leading to product B (cis-2-phenyl-2-butene) and 30 pathways leading to product C (3-phenyl-1-butene), respectively, with TSs located successfully. Interestingly, these products are in agreement with those from experimental data. No 19



pathways in which water molecule directly participated in the reaction have been located. To gain insight into the mechanism of this reaction, we have shown two typical reaction pathways in Figure 3, in which products B and C are formed respectively. These two pathways undergo a six-electron, five-membered cyclic transition structure, where a hydrogen atom transfer from carbon to oxygen is concerted with a bond cleavage of C-N bond. For pathway (1), product B was generated from the reactant (A1) via a concerted process (TSA1/B1) with the H13 transfer from C12 to O33 and the cleavage of bond C14-N24. The energy barrier for this step is 29.1 kcal/mol. In the pathway (2), a hydrogen atom (H23) transfer from the terminal methyl of the reactant A2 to O33 and the cleavage of bond C14-N24 occur simultaneously, leading to the product C1 with an energy barrier of 32.6 kcal/mol. By analyzing the structures of stationary points involved in these two pathways, we found that the oxygen atom of the amine oxide prefers to form two to three hydrogen bonds (O-HW) with the hydrogen atoms of surrounding water molecules. In addition, the hydrogen bonds in the reactant are shorter, and thus stronger than those in the corresponding TSs. Thus, the formed hydrogen bonds between the reactant and water molecules in the aqueous solution are unfavorable for this reaction.

20


Page 20 of 36

Page 21 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Two typical pathways leading to products cis-2-phenyl-2-butene (B, pathway (1)) and 3-phenyl-1-butene (C, pathway (2)) located by the MD/CD program for Cope elimination of amine oxide in aqueous solution. For simplicity, only nearest water molecules, which are forming hydrogen bond with the amine oxide, are illustrated among 100 H2O of the reactant/solvent cluster. Key distances (in Å) are labeled. Cartesian coordinates for all listed structures are available in SI.

It is of interest to compare the MD/CD results for the Cope elimination of amine oxide in aqueous and DMSO solution. The details of the MD/CD results in DMSO solution are available in the SI. A comparison of the zero-point energy corrected electronic energy barriers ΔE (kcal/mol) along dozens of pathways for this reaction in aqueous and DMSO solution is given in Figure 4. It shows that the energy barriers for 21



pathways in aqueous solution are, on average, higher than those in DMSO solution. The analysis of optimized structures for stationary points in aqueous and DMSO solution also shows that the reaction in DMSO solution is more favorable than that in aqueous solution. This qualitative conclusion is in good agreement with the experiment results and a previous theoretical study reported by Jorgensen group.63 Since the free energy barriers for this reaction in aqueous and DMSO solution have been reported previously,63 we do not carry out PMF calculations in this work.

Figure 4. Comparison of the electronic energy barriers ΔE (kcal/mol) corrected by zeropoint vibrational energy along dozens of pathways in the Cope elimination of amine oxide in aqueous solution (blue stars) and DMSO solution (red stars). 22


Page 22 of 36

Page 23 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.2 Dehydration of Methanediol in Aqueous Solution The solvent effect plays a very important role for the dehydration of methanediol in aqueous solution.64 Recent theoretical studies31,32 with a microsolvation model demonstrated that water molecules may act as a catalyst during the dehydration process of methanediol. Thus, the dehydration process of methanediol would be a good example to test the applicability of the refined MD/CD method. Here we first investigate the possible lower-barrier pathways for this reaction in aqueous solution with the refined MD/CD method ( Ecut

 50 kcal/mol). In addition to hydrogen atoms bonded to the carbon atom,

all other atoms of methanediol are chosen as target atoms. Then, we further perform PMF calculations along several reaction pathways to distinguish the most likely reaction pathway. Our MD/CD calculations have produced 104 reaction pathways leading to the dehydration product F, a formaldehyde and a water molecule, with TSs located successfully. According to the number of water molecules, which are directly involved in the reaction, these pathways can be classified into five groups (group 1 to 5 in Figure 5). Typical TS structures and pathways for each group are shown in Figure 5 and Figure S1, respectively. Along all pathways leading to product F, a hydrogen transfer and the cleavage of the C-O bond occur in a concerted way. For pathways in group 1 (no water molecules are involved), the product F was generated from a hydrogen transfer process (from one hydroxyl of the methanediol to another hydroxyl group) and a concomitant C-O bond 23



cleavage (via TS1). While for H2O-catalyzed pathways (in group 1 to 5), ring-like TS structures are formed between the methanediol and surrounding water molecules. These water molecules act as a bridge relaying the proton transfer from one hydroxyl of the methanediol to another hydroxyl. As shown in Figure 5, the electronic energy barriers in group 1 are on average much higher than those in other groups. These results indicate that the direct involvement of water molecules can reduce the energy barriers. Specifically, our results show that along pathways of group 3 (2H2O-catalyzed) and group 4 (3H2Ocatalyzed), the energy barriers are somewhat lower than those in other groups. Moreover, our calculations have produced 15 and 16 pathways in group 3 and 4, respectively, but only 3 pathways in group 2 (1H2O-catalyzed) and 2 pathways in group 5 (4H2O-catalyzed). It is worth noting that the conformers of the reactant are extracted from SITS-MD simulation rather than randomly selected, and all pathways are located by the MD/CD method automatically (without manually assigning the number of catalyzed water molecules). Thus, our calculations predict that the dehydration of methanediol may be catalyzed by two or three water molecules, as shown in pathways of group 3 and group 4. For these two possible reaction pathways located by the MD/CD method, we further calculate the corresponding free energy changes from PMF calculations. Besides the 2H2O-catalyzed and 3H2Ocatalyzed reaction modes, the reaction mode catalyzed by one water molecule have also been studied for comparison.

24


Page 24 of 36

Page 25 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 5. The zero-point energy corrected electronic energy barriers ΔE (kcal/mol) as a function of the number of catalyzed water molecules that are directly involved in the dehydration reaction of methanediol. Typical TSs for five groups are illustrated. For simplicity, only the catalyzed water molecules and other water molecules directly hydrogen bonded to methanediol in the reactant/solvent cluster are shown. Cartesian coordinates for all these typical TSs are available in SI.

The free energy barriers determined by PMF calculations for three reaction modes of the dehydration of methanediol in aqueous solution are summarized in Table 1. One can see that the 2H2O-catalyzed reaction mode has the lowest free energy barrier, 22.6±0.3 kcal/mol. This value is in good agreement with the experimental value (20.3±1.3 kcal/mol) 25



estimated from kd  4.96 10  e 7

6705/T -1 64

s

with the transition state theory. The free

energy barriers of the 1H2O-catalyzed or 3H2O-catalyzed reaction modes are 1.8~1.9 kcal/mol higher than that in the 2H2O-catalyzed pathway. For each reaction mode, the determined critical structures and the corresponding 2-D free energy surface are shown in Figure 6, Figure S3 and Figure S4. Thus, one can conclude that the most likely reaction mode for methonediol in aqueous solution is the dehydration of methonediol catalyzed by two water molecules.

Table 1. Free energy barriers, ΔG (kcal/mol, at 25C°), for three reaction modes of the dehydration of methanediol in aqueous solution determined by PMF calculations. In three reaction modes, the number of catalyzed water molecules is one, two, and three, respectively. Number of water molecules

ΔG

Experiment

1H2O

24.4±0.1

2H2O

22.6±0.4

3H2O

24.5±0.2

26


20.3 ±1.3

Page 26 of 36

Page 27 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6. The dehydration of methanediol in aqueous solution catalyzed by two water molecules. (a) Critical structures (in Å); (b) Two-dimensional free energy surface along two chosen reaction coordinates: the breaking C1−O6 distance (RC1=d(C1−O6)) and the difference between the breaking O4−H5 and forming O6−H10 distances (RC2= d(O4−H5)d(O6−H10)); (c) Computed O4(methanediol)…H (water) (O4-Hw) RDFs for the reactant (blue line) and TS (red line); (d) Computed H7(methanediol)-…O(water) (H7-Ow) RDFs for the reactant (blue line) and TS (red line). 27



In the dehydration of methanediol, another role of water molecules is to form hydrogen bond networks with the reactant and transition structures, which may have a significant effect on the reaction barrier. Here the RDFs, g (r ) , for the hydrogen bonds O4…Hw between the O4 atom of methanediol and hydrogen atoms of waters, and the hydrogen bonds H7…Ow between the H7 atom of methanediol and the oxygen atom of waters, are analyzed. For the 2H2O-catalyzed reaction mode (shown in Figure 6(c) and 6(d)), we found that for the O4…Hw hydrogen bonds the centers of the first peak of the corresponding g (r ) for the reactant and the transition structures are both at 1.75 Å. Integration of the first peaks to the minima near 2.3 Å reveals that the average number of O4…Hw hydrogen bonds is 0.90 (in the reactant) and 1.8 (in the transition structure). Similarly, our calculations show that the average number of H7…Ow hydrogen bonds is 0.99 (in the reactant) and 1.00 (in the transition structure). These results indicate that the hydrogen-bond interaction between methanediol and its neighboring water molecules stabilizes the transition structure more than the reactant, thus reducing the barrier to reaction. Similar trends can also be found for the 1H2O-catalyzed and 3H2O-catalyzed reaction modes. The detailed analysis can be found in the SI.

4. CONCLUSION

We have developed a refined MD/CD method for automatically searching reaction 28


Page 28 of 36

Page 29 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


pathways of chemical reactions in solution. In this refined MD/CD method, the SITS-MD simulation was employed to sample conformers of reactants in a realistic solution environment. In order to treat the solvent effect accurately, dozens of the reactant/solvent clusters (including the reactants and hundreds of neighboring solvent molecules) from SITS-MD simulations were used for the subsequent CD calculations with the ONIOMtype QM/MM method. For QM/MM calculations, the reaction center (consisting of the reactants and a few nearest-neighboring solvent molecules) is calculated at the QM level and the rest of the cluster is treated molecular mechanically. With this strategy, the MD/CD method is able to search reaction pathways, in which solvent molecules may directly participate. The performance of this refined MD/CD method was tested for two reactions: Cope elimination of threo-N, N-Dimethyl-3-phenyl-2-butylamine oxide in aqueous and DMSO solutions, and dehydration of methanediol in aqueous solution. For the Cope elimination reaction, our calculations can successfully locate pathways leading to the observed products, and the obtained energy barrier for this reaction in aqueous solution is found to be noticeably higher than that in DMSO. For the second reaction, our calculations can efficiently locate several reaction pathways, in which water molecules may act as catalysts. Based on the possible reaction modes from the MD/CD study, we have calculated PMF along different reaction coordinates to obtain the corresponding free energy barriers. Our calculations show that the most likely reaction mode for the second reaction corresponds 29



to the dehydration of methonediol catalyzed by two water molecules with a free energy barrier of 22.6±0.3 kcal/mol, which is in good agreement with the experimental value of 20.3±1.3 kcal/mol. Thus, the present study demonstrates that the refined MD/CD method is a promising tool for predicting low-energy pathways for reactions in solution. More work should be done in extending this method to more complicated reactions in different condensed phase environments.

ASSOCIATED CONTENT

Supporting Information The details of PMF calculations, scheme of geometrical comparison between two structures, typical pathways located by MD/CD method, cartesian coordinates for key structures, and PMF calculation results for dehydration of Methanediol in aqueous solution. This material is available free of charge via the Internet at http://pubs.acs.org.

AUTHOR INFORMATION

Corresponding Author *E-mail: [email protected]. Notes 30


Page 30 of 36

Page 31 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The authors declare no competing financial interest.

ACKNOWLEDGMENTS

The authors thank Prof. Jing Ma (Nanjing University) for her fruitful discussions. Manyi Yang was also indebted to Mr. Jingxiang Zou and Mr. Yuqi Wang for valuable discussions. This work was supported by the National Natural Science Foundation of China (Grant Nos. 21333004, 21673110, U1430237, 21573006). All calculations in this work have been done on the IBM Blade cluster system in the High Performance Computing Center of Nanjing University.

REFERENCES

(1) Houk, K. N.; Cheong, P. H., Computational Prediction of Small-Molecule Catalysts, Nature 2008, 455, 309-313. (2) Wang, G.; Zhang, H.; Zhao, J.; Li, W.; Cao, J.; Zhu, C.; Li, S., Homolytic Cleavage of a B-B Bond by the Cooperative Catalysis of Two Lewis Bases: Computational Design and Experimental Verification, Angew. Chem., Int. Ed. 2016, 55, 5985-5989. (3) Maeda, S.; Ohno, K.; Morokuma, K., Systematic Exploration of the Mechanism of Chemical Reactions: The Global Reaction Route Mapping (GRRM) Strategy Using the ADDF and AFIR Methods, Phys. Chem. Chem. Phys. 2013, 15, 3683-3701. (4) Sameera, W. M.; Maeda, S.; Morokuma, K., Computational Catalysis Using the Artificial Force Induced Reaction Method, Acc. Chem. Res 2016, 49, 763-773. (5) Maeda, S.; Ohno, K., Global Mapping of Equilibrium and Transition Structures on Potential Energy Surfaces by the Scaled Hypersphere Search Method: Applications to ab initio Surfaces of Formaldehyde and Propyne Molecules, J. Phys. Chem. A 2005, 109, 5742-5753. (6) Maeda, S.; Ohno, K., A New Method for Constructing Multidimensional Potential Energy Surfaces by a Polar Coordinate Interpolation Technique, Chem. Phys. Lett. 2003, 381, 177-186. (7) Ohno, K.; Maeda, S., Global Reaction Route Mapping on Potential Energy Surfaces of Formaldehyde, Formic Acid, and Their Metal-Substituted Analogues, J. Phys. Chem. A 2006, 110, 893331



8941. (8) Maeda, S.; Morokuma, K., Communications: A Systematic Method for Locating Transition Structures of A+B→X Type Reactions, J. Chem. Phys. 2010, 132, 241102. (9) Maeda, S.; Morokuma, K., Finding Reaction Pathways of Type A + B→X: Toward Systematic Prediction of Reaction Mechanisms, J. Chem. Theory Comput. 2011, 7, 2335-2345. (10) Maeda, S.; Harabuchi, Y.; Takagi, M.; Taketsugu, T.; Morokuma, K., Artificial Force Induced Reaction (AFIR) Method for Exploring Quantum Chemical Potential Energy Surfaces, Chem. Rec. 2016, 16, 2232-2248. (11) Maeda, S.; Taketsugu, T.; Morokuma, K., Exploring Transition State Structures for Intramolecular Pathways by the Artificial Force Induced Reaction Method, J. Comput. Chem. 2014, 35, 166-173. (12) Behn, A.; Zimmerman, P. M.; Bell, A. T.; Head-Gordon, M., Incorporating Linear Synchronous Transit Interpolation into the Growing String Method: Algorithm and Applications, J. Chem. Theory Comput. 2011, 7, 4019-4025. (13) Zimmerman, P. M., Automated Discovery of Chemically Reasonable Elementary Reaction Steps, J. Comput. Chem. 2013, 34, 1385-1392. (14) Shang, C.; Liu, Z. P., Stochastic Surface Walking Method for Structure Prediction and Pathway Searching, J. Chem. Theory Comput. 2013, 9, 1838-1845. (15) Zhang, X.-J.; Shang, C.; Liu, Z.-P., Double-Ended Surface Walking Method for Pathway Building and Transition State Location of Complex Reactions, J. Chem. Theory Comput. 2013, 9, 5745-5753. (16) Zhang, X. J.; Liu, Z. P., Reaction Sampling and Reactivity Prediction Using the Stochastic Surface Walking Method, Phys. Chem. Chem. Phys. 2015, 17, 2757-2769. (17) Martínez-Núñez, E., An Automated Method to Find Transition States Using Chemical Dynamics Simulations, J. Comput. Chem. 2014, 36, 222-234. (18) Martínez-Núñez, E., An Automated Transition State Search Using Classical Trajectories Initialized at Multiple Minima, Phys. Chem. Chem. Phys. 2015, 17, 14912-14921. (19) Cerjan, C. J.; Miller, W. H., On Finding Transition States, J. Chem. Phys. 1981, 75, 2800. (20) Basilevsky, M. V., The Topography of Potental-Energy Surfances, Chem. Phys. 1982, 67, 337346. (21) Koča, J., A Mathematical Model of the Logical Structure of Chemistry. A Bridge between Theoretical and Experimental Chemistry and a General Tool for Computer-Assisted Molecular Design, Theor. Chim. Acta. 1991, 80, 51-62. (22) Sun, J.-Q.; Ruedenberg, K., Gradient Extremals and Steepest Descent Lines on Potential Energy Surfaces, J. Chem. Phys. 1993, 98, 9707. (23) Fadrna.´, E.; Koča, J., Single-Coordinate-Driving Method Coupled with Simulated Annealing. An Efficient Tool to Search Cinformational Space, J. Phys. Chem. B 1997, 101, 7863-7868. (24) Cernohorsky, M.; Kettou, S.; Koca, J., Vader: New Software for Exploring Interconversions on Potential Energy Surfaces, J. Chem. Inf. Comput. Sci. 1999, 39, 705-712. (25) Kim, Y.; Choi, S.; Kim, W. Y., Efficient Basin-Hopping Sampling of Reaction Intermediates through Molecular Fragmentation and Graph Theory, J. Chem. Theory Comput. 2014, 10, 2419-2426. 32


Page 32 of 36

Page 33 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(26) Yang, M.; Zou, J.; Wang, G.; Li, S., Automatic Reaction Pathway Search Via Combined Molecular Dynamics and Coordinate Driving Method, J. Phys. Chem. A 2017, 121, 1351-1361. (27) Tomasi, J.; Persico, M., Molecular Interactions in Solution: An Overview of Methods Based on Continuous Distributions of the Solvent, Chem. Rev. 1994, 94, 2027-2094. (28) Tomasi, J.; Mennucci, B.; Cammi, R., Quantum Mechanical Continuum Solvation Models, Chem. Rev. 2005, 105, 2999-3094. (29) Fukaya, H.; Morokuma, K., A Theoretical Study of the Mechanism of Selective Fluorination of Saturated Hydrocarbons by Molecular Fluorine. Participation of CHCl3 Solvent Molecules in the Ionic Process, J. Org. Chem. 2003, 68, 8170-8178. (30) Goldsmith, B. R.; Hwang, T.; Seritan, S.; Peters, B.; Scott, S. L., Rate-Enhancing Roles of Water Molecules in Methyltrioxorhenium-Catalyzed Olefin Epoxidation by Hydrogen Peroxide, J. Am. Chem. Soc. 2015, 137, 9604-9616. (31) Inaba, S., Theoretical Study of Decomposition of Methanediol in Aqueous Solution, J. Phys. Chem. A 2015, 119, 5816-5825. (32) Inaba, S.; Sameera, W. M., Dehydration of Methanediol in Aqueous Solution: An Oniom(QM/MM) Study, J. Phys. Chem. A 2016, 120, 6670-6676. (33) Maeda, S.; Ohno, K.; Morokuma, K., An Automated and Systematic Transition Structure Explorer in Large Flexible Molecular Systems Based on Combined Global Reaction Route Mapping and Microiteration Methods, J. Chem. Theory Comput. 2009, 5, 2734-2743. (34) Maeda, S.; Abe, E.; Hatanaka, M.; Taketsugu, T.; Morokuma, K., Exploring Potential Energy Surfaces of Large Systems with Artificial Force Induced Reaction Method in Combination with ONIOM and Microiteration, J. Chem. Theory Comput. 2012, 8, 5058-5063. (35) Gao, Y. Q., An Integrate-over-Temperature Approach for Enhanced Sampling, J. Chem. Phys. 2008, 128, 064105. (36) Yang, L.; Qin Gao, Y., A Selective Integrated Tempering Method, J. Chem. Phys. 2009, 131, 214109. (37) Yang, L.; Liu, C. W.; Shao, Q.; Zhang, J.; Gao, Y. Q., From Thermodynamics to Kinetics: Enhanced Sampling of Rare Events, Acc. Chem. Res. 2015, 48, 947-955. (38) Zhang, J.; Yang, Y. I.; Yang, L.; Gao, Y. Q., Conformational Preadjustment in Aqueous Claisen Rearrangement Revealed by SITS-QM/MM MD Simulations, J. Phys. Chem. B 2015, 119, 5518-5530. (39) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.et al. Gaussian 09,Revision B.01; Gaussian, Inc: Wallingford, CT, 2009. (40) Tao, P.; Schlegel, H. B., A Toolkit to Assist Oniom Calculations, J. Comput. Chem. 2010, 23632369. (41) Chung, L. W.; Sameera, W. M.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.et al, The Oniom Method and Its Applications, Chem. Rev. 2015, 115, 5678-5796. (42) Bao, J. L.; Truhlar, D. G., Variational Transition State Theory: Theoretical Framework and Recent Developments, Chem. Soc. Rev. 2017, 46, 7548-7596. (43) Nieves-Quinones, Y.; Singleton, D. A., Dynamics and the Regiochemistry of Nitration of Toluene, J. Am. Chem. Soc. 2016, 138, 15167-15176. 33



(44) Wang, J. M.; 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) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L., Comparison of Simple Potential Functions for Simulating Liquid Water, J. Chem. Phys. 1983, 79, 926-935. (46) Case, D. A.; Betz, R. M.; Botello-Smith, W.; Cerutti, D. S.; Cheatham, T. E.; Darden, I., T.A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.et al. Amber 2016; University of California: San Francisco, 2016. (47) Bayly, C. I.; Cieplak, P.; Cornell, W.; Kollman, P. A., A Well-Behaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: The Resp Model, J. Phys. Chem. 1993, 97, 10269-10280. (48) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Kollmann, P. A., Application of Resp Charges to Calculate Conformational Energies, Hydrogen Bond Energies, and Free Energies of Solvation, J. Am. Chem. Soc. 1993, 115, 9620-9631. (49) Friesner, R. A.; Guallar, V., ab initio Quantum Chemical and Mixed Quantum Mechanics/Molecular Mechanics (QM/MM) Methods for Studying Enzymatic Catalysis, Annu. Rev. Phys. Chem. 2004, 56, 389-427. (50) Zhang, Y., Pseudobond ab initio QM/MM Approach and Its Applications to Enzyme Reactions, Theor. Chem. Acc. 2006, 116, 43-50. (51) Hu, H.; Yang, W., Free Energies of Chemical Reactions in Solution and in Enzymes with ab initio Quantum Mechanics/Molecular Mechanics Methods, Annu. Rev. Phys. Chem. 2008, 59, 573-601. (52) Gao, J.; Truhlar, D. G., Quantum Mechanical Methods for Enzyme Kinetics, Annu. Rev. Phys. Chem. 2002, 53, 467-505. (53) Patey, G. N.; Valleau, J. P., A Monte Carlo Method for Obtaining the Interionic Potential of Mean Force in Ionic Solution, J. Chem. Phys. 1975, 63, 2334-2339. (54) Boczko, E. M.; Brooks, C. L., Constant-Temperature Free Energy Surfaces for Physical and Chemical Processes, J. Phys. Chem. 1993, 97, 4509-4513. (55) Roux, B., The Calculation of the Potential of Mean Force Using Computer Simulations, Comput. Phys. Commun. 1995, 91, 275-282. (56) Shao, Y.; Molnar, L. F.; Jung, Y.; Kussmann, J.; Ochsenfeld, C.; Brown, S. T.; Gilbert, A. T.; Slipchenko, L. V.; Levchenko, S. V.; O'Neill, D. P.et al, Advances in Methods and Algorithms in a Modern Quantum Chemistry Program Package, Phys. Chem. Chem. Phys. 2006, 8, 3172-3191. (57) Case, D. A.; Betz, R. M.; Botello-Smith, W.; Cerutti, D. S.; Cheatham, T. E.; Darden, I., T.A.; Duke, R. E.; Giese, T. J.; Gohlke, H.; Goetz, A. W.et al. Amber 2012; University of California: San Francisco, 2012. (58) Chai, J. D.; Head-Gordon, M., Long-Range Corrected Hybrid Density Functionals with Damped Atom-Atom Dispersion Corrections, Phys. Chem. Chem. Phys. 2008, 10, 6615-6620. (59) Stephens, P. J.; Devlin, F. J.; Chabalowski, C. F.; Frisch, M. J., Ab initio Calculation of Vibrational Absorption and Circular Dichroism Spectra Using Density Functional Force Fields, J. Phys. Chem. 1994, 98, 11623-11627. (60) Cram, D. J.; McCarty, J. E., Studies in Sterochemistry. Xxiv. The Preparation and Determination of Configuration of the Isomers of 2-Amino-3-Phenylbutane, and the Steric Course of the Amine Oxide 34


Page 34 of 36

Page 35 of 36 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Pyrolysis Reaction in This System, J. Am. Chem. Soc. 1954, 76, 5740-5745. (61) Cram, D. J.; Sahyun, M. R. V., Room Temperature Wolff-Kishner Reduction and Cope Elimination Reactions, J. Am. Chem. Soc. 1962, 84, 1734-1735. (62) Sahyun, M. R. V.; Cram, D. J., Studies in Stereochemistry. Xxxv. Mechanism of Ei Reaction of Amine Oxides, J. Am. Chem. Soc. 1963, 85, 1263-1268. (63) Acevedo, O.; Jorgensen, W. L., Cope Elimination: Elucidation of Solvent Effects from QM/MM Simulations, J. Am. Chem. Soc. 2006, 128, 6141-6146. (64) Winkelman, J. G. M.; Ottens, M.; Beenackers, A. A. C. M., The Kinetics of the Dehydration of Methylene Glycol, Chem. Eng. Sci. 2000, 55, 2065-2071.

35



For Table of Contents use only

Combined Molecular Dynamics and Coordinate Driving Method for Automatic Reaction Pathway Search of Reactions in Solution Manyi Yang, † Lijiang Yang,‡ Guoqiang Wang,† Yanzi Zhou,*,† Daiqian Xie† and Shuhua Li*,† †School

of Chemistry and Chemical Engineering, Key Laboratory of Mesoscopic

Chemistry of Ministry of Education, Institute of Theoretical and Computational Chemistry, Nanjing University, Nanjing, 210023, People’s Republic of China ‡

Beijing National Laboratory for Molecular Sciences, College of Chemistry and

Molecular Engineering, Peking University, Beijing 100871, People’s Republic of China

36


Page 36 of 36

Combined Molecular Dynamics and Coordinate Driving Method for

Recommend Documents