Variational Flooding Study of a SN2 Reaction - The Journal of

Subscriber access provided by UNIV OF CALIFORNIA SAN DIEGO LIBRARIES

Letter N

Variational Flooding Study of a S 2 Reaction GiovanniMaria Piccini, James McCarty, Omar Valsson, and Michele Parrinello J. Phys. Chem. Lett., Just Accepted Manuscript • DOI: 10.1021/acs.jpclett.6b02852 • Publication Date (Web): 10 Jan 2017 Downloaded from http://pubs.acs.org on January 11, 2017

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 free 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 accessible to all readers and 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.

The Journal of Physical Chemistry Letters 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 13

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 Journal of Physical Chemistry Letters

Variational Flooding Study of a SN 2 Reaction GiovanniMaria Piccini,†,‡ James J. McCarty,†,‡ Omar Valsson,†,‡ and Michele Parrinello∗,†,‡ †Department of Chemistry and Applied Biosciences, ETH Zurich, c/o USI Campus, Via Giuseppe Buffi 13, CH-6900, Lugano, Ticino, Switzerland ‡Facoltá di Informatica, Instituto di Scienze Computationali, Universit’a della Svizzera italiana (USI), Via Giuseppe Buffi 13, CH-6900, Lugano, Ticino, Switzerland E-mail: [email protected]

1

ACS Paragon Plus Environment


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 ACS Paragon Plus Environment

Page 2 of 13

Page 3 of 13

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 study of reaction dynamics is key to understanding the nature of chemical transformations. The textbook approach to study chemical reactions starts from computing the potential energy surface (PES) on which reactants and products are made to move either using classical or quantum mechanics. 1 Alternatively, theories like transition state (TS) are used to compute rates. A precise determination of the rates necessitates the use of a high level theory, furthermore an accurate representation of the PES becomes more and more challenging as the number of degrees of freedom increases. Although this approach leads to very accurate results, it is computationally expensive and only relatively small systems can be studied in this way. For this reason it has been suggested to calculate the forces on the fly. 2 In doing so, often the accuracy of the interaction forces is reduced due to computational restrictions. This loss of accuracy is compensated by the fact that rather more complex and realistic systems can be studied. Ideally one would like to simulate the spontaneous evolution of the system from educts to products under the action of the interatomic forces. Unfortunately in most cases educts and products are separated by large barriers that cannot be overcome in a direct simulation. Therefore, special methods are needed. If the transition state is known, a variety of methods allow the reaction path to be calculated. 3 However, in complex systems the TS is not easy to be determined and its notion needs to be generalized. 4 Based on previous works by Grubm¨ uller 5 and Voter 6 we have recently proposed a method to circumvent this problem and have combined their ideas with two enhanced sampling methods developed in our group, namely metadynamics 7,8 and variationally enhanced sampling (VES). 9 This has lead to two methods for calculating rates in a rare event scenario. We shall call these two methods infrequent metadynamics 10 and variational flooding. 11 Both approaches have proven to be useful in several applications and have allowed rates in complex systems to be computed. More recently the application of infrequent metadynamics to – −− the symmetric CH3 Cl + Cl – ↽ −⇀ − CH3 Cl + Cl SN 2 reaction has been reported, demon3



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

Page 4 of 13

strating its potential in the study of chemical reactions. 12 Here, using instead variational – flooding, we shall study a SN 2 reaction as well, namely CH3 F + Cl – − ↽− −⇀ − CH3 Cl + F . This reaction differs from that of Ref. 12 by having a highly asymmetric energy profile. The method allows to harvest a large number of reaction trajectories and study the reaction rate as a function of temperature. In this scheme all anharmonic effects are automatically included. Both metadynamics and VES are based on the introduction of a set of CVs s = s(R) where R are the atomic coordinates. We associate to the chosen set of CVs a free-energy surface (FES) 1 F (s) = log β

Z

dR δ (s − s(R)) e−βU (R)

(1)

where U (R) is the interaction potential, and β = (kB T )−1 is the inverse temperature. Here and in the following we shall drop immaterial constants. As in umbrella sampling a bias V (s(R)) that depends on the CV is introduced. The function of V (s) is to facilitate transitions from one metastable state to another. In VES the bias is constructed by minimizing the functional: 1 Ω [V ] = β

R

ds e−β[F (s)+V (s)] R + ds e−βF (s)

Z

ds p(s)V (s)

(2)

where p(s) is a chosen target probability distribution that is assumed to be normalized to one. It is easy to show that the bias potential that minimizes Ω [V ] is related to F (s) by the relation

V (s) = −F (s) −

1 log p(s) β

(3)

and that such a solution corresponds to the global minimum being Ω [V ] a convex functional. 9 The variational property of Ω [V ] is exploited by assuming a functional form of the bias potential V (s; α) that depends on a set of variational parameters α = (α1 , α2 , . . . , αK ). For 4


Page 5 of 13

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


instance, the bias potential can be linearly expanded as

V (s; α) =

X

αk fk (s)

(4)

k

where the fk (s) form an appropriate set of basis functions such as plane waves or Chebyshev polynomials, however, other choices for V (s; α) are possible. 13 Having defined the bias in this way the functional can be minimized numerically (see Ref. 9 for details) to obtain the corresponding FES (see Eq. 3). The idea behind variational flooding is the following. Having the variational principle for Ω [V ] we calculate a bias that depends on s that is able to facilitate the escape from a deep basin. In the construction of the potential a free-energy cut-off on V (s) is imposed such that the bias is zero in the region in which the system transits for one basin to another. Under these conditions, using the arguments of Grubm¨ uller 5 and Voter, 6 it can be shown that the escape times measures in the biased MD are related to the physical ones by

t∗ = tM D heβV (s) iV .

(5)

where t∗ is physical time, tM D is the simulation time in the biased simulation, and the average heβV (s) iV , obtained over the biased sample is the acceleration factor that measures how much the escape time is boosted due to the presence of the bias V (s). The validity of the rare event assumption requires that the exit times must have a Poisson distribution. Thus the cumulative distribution function (CDF) is fitted to the Poisson expression P = 1 − exp(t/τ ) and the statistical validity of this fit measured by the p-values test 14 . In this study we considered the asymmetric SN 2 nucleophilic substitution reaction of – fluoromethane and chloromethane: 15,16 CH3 F + Cl – − ↽− −⇀ − CH3 Cl + F . We shall refer to the educt as A and the product as B. Such a reaction presents a very high barrier when substituting the fluorine atom with the chlorine one, making it particularly challenging 17 . In order to drive the reaction we use as CV the variable s = d1 − d2 where d1 and d2 are

5



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 C−F and C−Cl distances respectively. This CV is a good indicator of whether or not the reaction has taken place and reflects the concerted mechanism of the reaction, but it is only an approximated representation of the true reaction coordinate. Nevertheless since the method is very forgiving in this respect such a simple CV suffices for our purposes. We stress that although the chosen CV is one dimensional the system is treated in all its complexity and all degrees of freedoms are taken fully into account. The simulations have been performed using the semi-empirical model PM6 18 to ensure a qualitatively correct description of the bond breaking and forming process allowing long simulation times at low computational cost. the simulations have been carried out using the CP2K program 19 supplemented by a, soon to be released, private version of the PLUMED 2 code. 20 The system has been put in a 10 × 10 × 10 ˚ A supercell without periodic boundary conditions converging the energy at each self-consistent field iteration with a threshold accuracy of 10−7 Hartree. The time step used is 0.5 fs. To control the temperature the systems has been coupled to the velocity rescaling thermostat of Bussi et al. 21 Here the thermostat aims at mimicking in a simple manner the effect of the environment. We studied the reaction at three different temperatures, namely 600, 900, and 1200 K. Harmonic restraining walls have been placed on d1 and d2 for distances larger than 3 ˚ A to prevent the halide anions from moving too far from the organic substrate. The bias potential is represented using an expansion in Chebyshev polynomials up to order 30 to obtain the FES profile and up to 50 to obtain the flooding potential. A cut-off of 133 kJmol−1 has been imposed on the height of the bias potential. During the simulation the bias is optimized by varying the coefficients of the expansion every 1 ps. The step size used in the average stochastic gradient descent 22 was first set to 10 kJmol−1 to allow a full exploration of the CV space. Subsequently the simulations have been restarted reducing such step size to 1 kJmol−1 in order to refine the coefficients. To study the transitions between the reactants and products, 30 independent trajectories have been evolved all initiated after thermalization in the A basin. During the simulation the

6


Page 6 of 13

Page 7 of 13

system visits several times the two basins. These long trajectories have been subsequently separated into segments each describing one single reactive event either A to B or B to A. To ensure that the system has had time to thermalize after the transition in the new basin, trajectory segments shorter than 1 ps have been discarded. We first calculated the free-energy associated to s (see Figure 1). It shows a deep minimum for the educt and a much shallower one for the product. This FES structure suggests the use of a free-energy cut-off that is lower then the B basin minimum, yet sufficiently high for the escape events from the biased basin to be frequent. The choice of a cut-off of 133 kJmol−1 ensures this to be the case. In such a manner we could continue the simulation for a very long time and observe hundreds of transitions, thus reducing the statistical error considerably. 180 A

160

B

140 Energy (kJ mol−1 )

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


120 100 80 60 40 U (s) F (s) F (s) + V (s)

20 0 -3

-2

-1

0 s (˚ A)

1

2

3

Figure 1: Energy profiles along the CV. Potential energy surface, U (s), free-energy surface, F (s), and flooded FES, F (s) + V (s), all at 1200 K. Note at this temperature the sizable entropy effects. In this way we impose a dynamical asymmetry. the transitions A → B are stimulated by the presence of the bias while the reverse B → A transitions occur spontaneously. This implies that the measured exit times A → B need to be renormalized according to Eq. 5 7



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

Page 8 of 13

while the B → A times have their physical value. In order to apply the theory on which Eq. 5 is based 11 it is necessary to verify that even after the imposition of the bias one is still in the rare event regime. If this condition is satisfied, the distribution of escape times must be Poissonian and the statistical accuracy of this assumption verified 14 . This is shown in Figure 2 and Table 1. It is seen that indeed the A → B transitions have the desired behaviour. In the B → A transition small deviations from a perfect Poisson behaviour are observed at short times. However, these small deviations are of little practical consequence and the overall behaviour is very close to being Poissonian as reflected by the p-values (see Table 1). Table 2 reports the rate constants calculated using the Poisson fit procedure for both reaction directions. Both kA→B and kB→A calculated in different ways exhibit the expected Arrhenius behavior (see Figure 3). From the Arrhenius plot the activation energies can be calculated (see Table 3). these values compare well with the potential energy barrier. Table 1: Transition times, τ and Standard Error (all in s−1 ), total number of transition observed, n, and p-value for both reactions at 600, 900, and 1200 K. The latter quantity has been calculated according to Ref. 14 for a randomly selected subset of 100 escape times at each temperature. The p-test, assessing the reliability of the reconstructed dynamics, resulted to be always positive (p > 0.05) for both reaction directions at the three different temperatures. τ

Standard Error

n

p

A→B

600 K 900 K 1200 K

3.18 × 101 1.08 × 10−3 6.48 × 10−6

1.06 2.49 × 10−5 1.14 × 10−7

814 2045 1259

0.35 0.83 0.39

B→A

600 K 900 K 1200 K

1.72 × 10−11 9.43 × 10−12 5.64 × 10−12

6.97 × 10−13 1.87 × 10−13 1.63 × 10−13

897 2390 1309

0.84 0.15 0.65

8


Page 9 of 13

1 Cumulative Probability

A→B 0.8 0.6 −→ 0.4 0.2 0

10−15

101

10−7 t∗ (s)

1 B→A Cumulative Probability

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


0.8 0.6 0.4 0.2 0

10−13

10−12 t∗ (s)

10−11

10−10

Figure 2: CDF and relative fit to a Poisson distribution function of the transition times for both reaction directions at 600 (green), 900 (blue), and 1200 K (purple). In the upper panel, the dashed lines are the distributions corresponding to the unscaled escape times.

9



Table 2: Rate constants, k, for both reaction directions at 600, 900, and 1200 K calculated using the Poisson fit of the escape times and the integral of the population time correlation function. All in s−1 .

600 K 900 K 1200 K

kA→B

kB→A

3.14 × 10−2 9.26 × 102 1.54 × 105

5.81 × 1010 1.06 × 1011 1.77 × 1011

ln(kA→B )

13.0

7.0

1.0

−5.0

A→B

26.0 ln(kB→A )

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

25.5

25.0 B→A 1.1 × 10−4 1.5 × 10−4 1.9 × 10−4 1/RT (kJ−1 mol) Figure 3: Arrhenius plots for both the reaction directions calculated using different methods.

10


Page 10 of 13

Page 11 of 13

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


Table 3: Activation energies calculated from Arrhenius plot fitting procedure, EAf it , and from PM6 PES, EAP M 6 , for both the reaction directions. All in kJmol−1 .

A→B B→A

EAf it

EAP M 6

153.75 11.90

157.22 9.32

Variationally enhanced sampling is a powerful tool to explore free-energy landscapes characterized by long-lived metastable states. The general formalism of the method allows an ad hoc design of a free-energy flooding bias potential used to calculate transition times and, therefore, rate constants. In this work we have shown that the method can be applied to challenging cases such as chemical reactions for which the barrier is both high and steep. We conclude that the method is a promising tool for studying chemical reactions and can be applied to a large variety of chemical problems from surface science to biological systems.

Acknowledgement This research was supported by the European Union Grant No. ERC-2014-AdG-670227/VARMET and by the NCCR MARVEL, funded by the Swiss National Science Foundation. Calculations were carried out on the Mönch cluster at the Swiss National Supercomputing Center (CSCS).

References (1) Xie, J.; Hase, W. L. Rethinking the SN2 reaction. Science 2016, 352, 32–33. (2) Car, R.; Parrinello, M. Unified Approach for Molecular Dynamics and DensityFunctional Theory. Phys. Rev. Lett. 1985, 55, 2471–2474. (3) Truhlar, D. G.; Steckler, R.; Gordon, M. S. Potential energy surfaces for polyatomic reaction dynamics. Chem. Rev. 1987, 87, 217–236. 11



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

(4) Chandler, D. Statistical mechanics of isomerization dynamics in liquids and the transition state approximation. J. Chem. Phys. 1978, 68, 2959–2970. (5) Grubm¨ uller, H. Predicting slow structural transitions in macromolecular systems: Conformational flooding. Phys. Rev. E 1995, 52, 2893–2906. (6) Voter, A. F. Hyperdynamics: Accelerated Molecular Dynamics of Infrequent Events. Phys. Rev. Lett. 1997, 78, 3908–3911. (7) Laio, A.; Parrinello, M. Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562–12566. (8) Valsson, O.; Tiwary, P.; Parrinello, M. Enhancing important fluctuations: Rare events and metadynamics from a conceptual viewpoint. Ann. Rev. Phys. Chem. 2016, 67, 159–184. (9) Valsson, O.; Parrinello, M. Variational Approach to Enhanced Sampling and Free Energy Calculations. Phys. Rev. Lett. 2014, 113, 090601. (10) Tiwary, P.; Parrinello, M. From Metadynamics to Dynamics. Phys. Rev. Lett. 2013, 111, 230602. (11) McCarty, J.; Valsson, O.; Tiwary, P.; Parrinello, M. Variationally Optimized FreeEnergy Flooding for Rate Calculation. Phys. Rev. Lett. 2015, 115, 070601. (12) Fleming, K. L.; Tiwary, P.; Pfaendtner, J. New Approach for Investigating Reaction Dynamics and Rates with Ab Initio Calculations. J. Phys. Chem. A 2016, 120, 299– 305, PMID: 26690335. (13) McCarty, J.; Valsson, O.; Parrinello, M. Bespoke Bias for Obtaining Free Energy Differences within Variationally Enhanced Sampling. J. Chem. Theory and Comput. 2016, 12, 2162–2169.

12


Page 12 of 13

Page 13 of 13

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


(14) Salvalaglio, M.; Tiwary, P.; Parrinello, M. Assessing the Reliability of the Dynamics Reconstructed from Metadynamics. J. Chem. Theory Comput. 2014, 10, 1420–1425, PMID: 26580360. (15) Szabó, István and Császár, Attila G. and Czakó, Gábor, Dynamics of the F- + CH3Cl [rightward arrow] Cl- + CH3F SN2 reaction on a chemically accurate potential energy surface. Chem. Sci. 2013, 4, 4362–4370. (16) Szabó, I.; Czakó, G. Revealing a double-inversion mechanism for the F-+ CH3Cl SN2 reaction. Nat. Commun. 2015, 6 . (17) Mugnai, M.; Cardini, G.; Schettino, V. An ab initio molecular dynamics study of the SN2 reaction F+CH3Cl → CH3F+Cl. J. Chem. Phys. 2003, 118, 2767–2774. (18) Stewart, J. J. Optimization of parameters for semiempirical methods V: modification of NDDO approximations and application to 70 elements. J. Mol. Model. 2007, 13, 1173–1213. (19) Hutter, J.; Iannuzzi, M.; Schiffmann, F.; VandeVondele, J. cp2k: atomistic simulations of condensed matter systems. WIREs Comput Mol Sci 2014, 4, 15–25. (20) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. Comput. Phys. Commun. 2014, 185, 604–613. (21) Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J. Chem. Phys. 2007, 126, 014101. (22) Bach, F.; Moulines, E. In Advances in Neural Information Processing Systems 26 ; Burges, C., Bottou, L., Welling, M., Ghahramani, Z., Weinberger, K., Eds.; Curran Associates, Inc., Red Hook, NY, 2013; pp 773–781.

13


Variational Flooding Study of a SN2 Reaction - The Journal of

Recommend Documents