Modeling Matter: Calculating Reaction Rates & Defining Coordinates in Complex Systems, Study notes of Chemistry

Download Modeling Matter: Calculating Reaction Rates & Defining Coordinates in Complex Systems and more Study notes Chemistry in PDF only on Docsity! NC STATE UNIVERSITY CHE596M Multi-Scale Modeling of Matter Dr. Erik E. Santiso Lecture 26: Chemical Reactions and Activated Processes (II) NC STATE UNIVERSITY Outlook • Finding reaction mechanisms: – The Zero-Temperature String Method • Calculating reaction rates and defining reaction coordinates in complex systems: – Introduction – The Bennett-Chandler method – “Blue Moon” molecular dynamics – Transition path sampling – Aimless Shooting + Likelihood Maximization – Finite Temperature String Method – Other methods NC STATE UNIVERSITY • The Zero-Temperature String Method [15, 16]: – It can be shown that NEB is a particular case of the zero- temperature string method. However, it is simpler and more efficient to parameterize the reaction path using arc length (intrinsic reaction coordinate) or energy-weighted arc length. – Example: http://www.cims.nyu.edu/~eve2/string.html Finding reaction mechanisms NC STATE UNIVERSITY Outlook • Finding reaction mechanisms: – The Zero-Temperature String Method • Calculating reaction rates and defining reaction coordinates in complex systems: – Introduction – The Bennett-Chandler method – “Blue Moon” molecular dynamics – Finite Temperature String Method – Transition path sampling – Aimless Shooting + Likelihood Maximization – Other methods NC STATE UNIVERSITY Chemical Reaction • A process that results in the interconversion of chemical species. • A chemical species is defined as “An ensemble of chemically identical molecular entities that can explore the same set of molecular energy levels on the time scale of the experiment” (From IUPAC glossary of terms used in physical organic chemistry, available online at http://www.chem.qmul.ac.uk/iupac/gtpoc/) • Key point: Species are defined in terms of a time scale. So, which chemical species are you?I don’t know yet. Give me some time… Introduction: Remember this slide? NC STATE UNIVERSITY • Reaction rates in chemistry are defined in terms of rate laws. Rate laws are phenomenological equations and cannot be derived. They arise as the result of experimental observations. We can find statistical mechanical expressions for phenomenological coefficients by assuming that the rate law is correct and then relating it to microscopic variables. 0 0.5 1 1.5 20 1 2 Time (min) C on ce nt ra tio n (m ol /L ) Hmm…Looks like AAA Ckr −= Reaction rates NC STATE UNIVERSITY Calculation of rate constants • The rate law for a chemical reaction is a phenomenological equation. We need to assume that it is correct and find an expression for the rate constant. • The derivation is in many ways similar to the derivation of the fluctuation-dissipation theorem we saw for Langevin dynamics, or Einstein’s equation for the diffusion coefficient. They relate the (linear) response of the system to a perturbation to the equilibrium fluctuations of microscopic quantities over time. kT2 2σ =γ ( ) ( )[ ]20lim6 rtr dt dD t −= ∞→ • We want to find a relation of this kind between the rate constant and the equilibrium fluctuations of the concentration… NC STATE UNIVERSITY Calculation of rate constants • Let us consider a simple reaction, , for which we know that the rate equations are: BA↔ BrAf B BrAf A ckck dt dc ckck dt dc −= +−= (26.1) • Where kf and kr are the forward and backward rate constants, and cA, cB are the number densities of A and B. We will assume that they are normalized, i.e. • At equilibrium, the concentrations are constant, and we have: 1=+ BA cc r f eq A eq B eq k k c cK == (26.2) NC STATE UNIVERSITY • The reaction coordinate is a function of the configuration. The value corresponds to the transition state between A and B, and defines a “dividing surface” in configurational space that separates reactants from products. If we now average over the canonical ensemble, we obtain the probability that the configuration is a “reactant” configuration, which is equal to cA: Calculation of rate constants ξ ‡ξ Aθ ( ) ( ) A A A cQ Hrd = β−θΓ =θ ∫ exp (26.7) • Here H is the Hamiltonian of the system. Although we are deriving the expression for the canonical ensemble, the final result is actually the same for other ensembles as well. For derivations see [2,3]. NC STATE UNIVERSITY Calculation of rate constants • Now we would like to find an expression for , the relaxation of the concentration upon an initial perturbation, so that we can relate it to equation (30.4). If we perturb the Hamiltonian so that the states with have a lower energy, the equilibrium concentration of A will be higher. We can do this by lowering the energy of these states: ( )tcA∆ ‡ξ<ξ AHH λθ−= 0 (26.8) Unperturbed Hamiltonian Small perturbation • If the perturbation is small, we can expand on a Taylor series and keep only the leading term: Ac∆ ( ) ( )2 0 0 0 λ+⎟ ⎠ ⎞ ⎜ ⎝ ⎛ λ∂ ∆∂ λ+=λ∆=θ−θ=∆ =λ λ Occc AAAAA (26.9) 0= 0≈ NC STATE UNIVERSITY Calculation of rate constants • Thus we have: [ ] [ ] [ ] [ ] [ ][ ] [ ][ ] [ ] 2 00 2 0 2 22 0 00 0 0 0 0 0 0 AAH H A H H A H H AAA A A A A A A A ed ed ed ed ed edccc θ−θβλ= ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ Γ βθΓβ − Γ βθΓ λ= = ⎪⎭ ⎪ ⎬ ⎫ ⎪⎩ ⎪ ⎨ ⎧ Γ θΓ λ∂ ∂ λ=⎟ ⎠ ⎞ ⎜ ⎝ ⎛ λ∂ ∂ λ=⎟ ⎠ ⎞ ⎜ ⎝ ⎛ λ∂ ∆∂ λ=∆ =λ λθ−β− λθ−β− λθ−β− λθ−β− =λ λθ−β− λθ−β− =λ=λ ∫ ∫ ∫ ∫ ∫ ∫ (26.10) • Or, since : [ ] [ ] ( )01 00 2 00 A eq B eq AAAAAA cccc ∆=βλ=θ−θβλ=θ−θβλ=∆ (26.11) AA θ=θ 2 NC STATE UNIVERSITY Calculation of rate constants • Now, from equation (26.11), we had: ( ) ( )eq B eq A Aeq B eq AA cc cccc 0 0 ∆=βλ⇒βλ=∆ (26.14) • Substituting this into (26.13) we get: ( ) ( ) ( ) ( ) ( )eq B eq A eq AAA AA cc ct ctc 2 0 0 0 −θθ ∆=∆ • From the phenomenological rate law we had: ( ) ( ) [ ] f eq B RRAA k ctctc =ττ−∆=∆ , exp0 NC STATE UNIVERSITY Calculation of rate constants • Comparing the last two equations, we get: ( ) ( ) ( ) ( )eq B eq A eq AAA R cc ct t 2 0 0 exp −θθ =τ− (26.15) • This equation is actually valid only for long times. As is the case with Einstein’s equation for the diffusion coefficient, and other autocorrelation function expressions, the short-time behavior of the function can be quite different. • Differentiating (26.15) with respect to time, we get: ( ) ( ) ( ) eq B eq A AA R R cc t t 0 0 exp1 θθ =τ− τ − (26.16) NC STATE UNIVERSITY Calculation of rate constants • If we now assume that the reaction is truly a rare event, the time sampled by our simulation will be much smaller than . Then we can assume that and we get the very important result: Rτ ( ) ( ) ( ) ( ) eq B eq A AA eq B eq A AA R cc t cc t 00 001 θθ −= θθ = τ − ( ) 1exp ≈τ− Rt (26.17) • The last equality can be obtained by changing time t to and reversing time. This only changes the sign of the autocorrelation function. If we now recall that we get the Bennett- Chandler equation: 0=t f eq BR kc=τ ( ) ( ) eq A AA f c t k 0 0 θθ = (26.18)for Rmol t τ<<<<τ NC STATE UNIVERSITY Calculation of rate constants • Using this result, we can factor the Bennett-Chandler equation as: ( )[ ] ( ) ( )[ ] ( )[ ] ( ) ( )[ ] ( )[ ] ( ) ( )[ ] TSTeq A f kc t k κ= ξθξξ−ξδ ξθξξ−ξδ ξ−ξθξξ−ξδ = 000 000 00 ‡ ‡ ‡‡ (26.23) • Here is the transmission coefficient. This is analogous to equation (26.20) in that one factor is an equilibrium quantity and the other is a dynamic quantity. • The TST constant assumes that everything that reaches the top of the barrier ultimately reacts. The transmission coefficient accounts for the possibility of recrossings of the barrier, and is always less than or equal to 1. • You may have seen equation (26.23) with an extra factor, which accounts for quantum tunneling through the barrier. This is a non- classical term and cannot be obtained from this theory. κ NC STATE UNIVERSITY The Bennett-Chandler Method • There are two ways we can use the Bennett-Chandler equation to calculate rate constants. • If we use equation (26.20), we can obtain the second (equilibrium) term by realizing that: [ ] { } { }barrier theof side reactant"" theon is system y that theProbabilit barrier theof topat the is system that the(density)y Probabilit‡ = θ ξ−ξδ A • The calculation of probabilities is equivalent to the calculation of free energies (remember the coarse-grained Hamiltonian?). There are standard techniques to deal with this problem, such as umbrella sampling or thermodynamic integration. We will see shortly an example of how to do this…. NC STATE UNIVERSITY The Bennett-Chandler Method • The first term in equation (26.20) is a conditional average: ( )[ ] ( ) ( )[ ] ( )[ ]‡ ‡‡ 0 00 ξ−ξδ ξ−ξθξξ−ξδ t Count only trajectories that start at the top of the barrier • We can obtain this quantity by using molecular dynamics. We sample only dynamic trajectories that start at the top of the energy barrier. This condition, however, introduces a bias in the initial distribution function, and we will need to correct for this. NC STATE UNIVERSITY TST + RxMC: Example • Decomposition of HI in the bulk, in solution and in carbon nanopores [6] NC STATE UNIVERSITY Outlook • Finding reaction mechanisms: – The Zero-Temperature String Method • Calculating reaction rates and defining reaction coordinates in complex systems: – Introduction – The Bennett-Chandler method – “Blue Moon” molecular dynamics – Transition path sampling – Aimless Shooting + Likelihood Maximization – Finite Temperature String Method – Other methods NC STATE UNIVERSITY Blue Moon Molecular Dynamics • From the point of view of the molecules, reactions happen only “once in a blue moon”. (actually blue moons are comparatively much more common, they happen every two or three years, not 60000! Last one was May 30, 2007, the next one will be December 2009). (Picture from Michael Myers’ blue moon web page at http://www.netaxs.com/~mhmyers/blue/blue.html) • The Blue Moon MD method is a very elegant way to obtain both terms in equation (26.20) from molecular dynamics. The main idea is to: (1) Use constrained dynamics to generate configurations at the top of the barrier + thermodynamic integration to obtain the probability ratio , then (2) Use these configurations as starting points to run unconstrained dynamics and estimate the dynamical term [ ] Aθξ−ξδ ‡ ( )[ ] ( ) ( )[ ] ( )[ ]‡‡‡ 000 ξ−ξδξ−ξθξξ−ξδ t NC STATE UNIVERSITY Blue Moon Molecular Dynamics [ ] ( ) ( ) ( ) ∫ ∫ β− β− θ ξ−ξδ =ξξ= θ ξ−ξδ U U erd erd PP A N N A * *‡ ‡ h wit (26.26) • We can differentiate with respect to and change variables to a new coordinate system in which is one of the coordinates (e.g. the eigenbasis of the Hessian). This way we obtain: ( )*ln ξP *ξ ( ) ( ) ( ) ( ) ( ) C C Z J ZJ P 21 21 * * * * lnln ln − − ξ∂ β−∂ = ξ−ξδ ξ∂ β−∂ ξ−ξδ = ξ∂ ξ∂ UU (26.27) • Here J is the Jacobian of the coordinate transformation. This equation can be integrated numerically to get .( )‡ξP ξ • The “equilibrium” part of the rate constant can be obtained using thermodynamic integration. First we rewrite it as: NC STATE UNIVERSITY Blue Moon Molecular Dynamics • In summary, a Blue Moon MD simulation goes as follows: (1) Run constrained dynamics with values of ranging from (the reactant side) to (the transition state). Use these simulations to evaluate the averages on the right-hand side of eq. (26.27). Integrate this equation numerically to find the “equilibrium” part of the rate constant. (2) Using as initial conditions the configurations from the constrained MD with , run “normal” MD simulations to calculate the averages on the right-hand side of eq. (26.24). This gives the “dynamic” part of the rate constant. (3) Multiply the two terms and voila!, you have the rate constant. ξ Aξ ‡ξ ‡ξ=ξ NC STATE UNIVERSITY Blue Moon MD: Example • Simulation of the SN2 reaction using ab initio molecular dynamics in the Blue Moon ensemble [11] FCHClClCHF 33 +→+ −− NC STATE UNIVERSITY Transition Path Sampling • In many cases, we simply cannot define a reaction coordinate because the process is too complex (e.g. nucleation, protein folding). We need a method to evaluate rate constants that does not require introducing a reaction coordinate. • From the Bennett-Chandler method, we saw that the dynamic part of the rate constant can be obtained by sampling only over reactive trajectories. If we know what the “reactants” and “products” are, we only need trajectories connecting the two. (Figure from ref. [12]) NC STATE UNIVERSITY Transition Path Sampling • We need to write the Bennett-Chandler equation in a way that does not require defining a reaction coordinate. One possibility is: ( ) ( ) ( ) ( ) A BA f t tCtC dt dk θ θθ == 0 with (26.28) • The correlation function counts how many trajectories that start at the reactants at end up being products at time t. This still requires us to define what the reactants and products are, but not a reaction coordinate. The correlation function can be rewritten as a ratio of path integrals: ( )tC 0=t ( ) ( ) ( ) ( ) ( )[ ] ( ) ( ) ( )[ ]∫ ∫ θ θθ = trtr trttr tC N A N N BA N PD PD 0 0 (26.29) NC STATE UNIVERSITY Transition Path Sampling • The integration in (26.29) is over the space of all dynamic trajectories. This is a bit different from the sampling over configurations that we are used to. The Transition Path Sampling (TPS) method attempts to evaluate , and hence the rate constant, by doing Monte Carlo sampling on the space of trajectories that connect reactants with products. ( )tC • A trajectory is in principle an infinite collection of configurations. In practice, it is approximated by a discrete set of points in configuration space . The probability of a given trajectory, is calculated differently depending on the dynamic method use to generate it (e.g. Langevin dynamics, “standard” MD, kinetic Monte Carlo). It is usually written as a product of transition probabilities from one point in the trajectory to the next, times the probability density of the initial point. ( ) ( ) ( )NMMNN rttrtrt ,,,,,,0 1100 == ( )[ ]tr NP NC STATE UNIVERSITY Transition Path Sampling How do we sample trajectory space? • “Shooting” moves: Change the momenta a little bit at an intermediate point, and integrate forward and backwards in time. • There are several possible “Monte Carlo moves” we can use to generate new trajectories. The two most common are: (From [12]) NC STATE UNIVERSITY Transition Path Sampling How do we sample trajectory space? • There are several possible “Monte Carlo moves” we can use to generate new trajectories. The two most common are: • “Shifting” moves: Delete a part of the trajectory at the beginning (or end) and integrate forward (backwards) in time from the end (beginning). (From [12]) NC STATE UNIVERSITY • There are other possible trial moves, e.g. configurational bias-like moves, to generate new trajectories. So far the combination of shooting and shifting seems to be the most efficient. • Note that we did not get away with calculating the reaction rate without finding free energies. In the equation for the reaction rate in the transition path ensemble, both the numerator and denominator are probabilities: Transition Path Sampling ( ) ( ) ( ) ( ) ( )[ ] ( ) ( ) ( )[ ] ( ) A AB N A N N BA N Z tZ trtr trttr tC = θ θθ = ∫ ∫ PD PD 0 0 • This can be related to the (average) reversible work needed to take a trajectory that ends anywhere, and make it end in the “products” region. This is in practice obtained using thermodynamic integration in trajectory space. NC STATE UNIVERSITY Aimless Shooting and Likelihood Maximization • As we discussed earlier, reaction coordinates are often found through trial-and-error. • The problem with the trial-and-error approach is that it requires either looking for hysteresis in a free energy profile, or generating a committor probability distribution for the trial reaction coordinate, both of which are very computationally expensive. • Alternatively, one can generate an unbiased set of reactive trajectories, and analyze the data afterwards to determine suitable reaction coordinates. This is the idea of the Aimless Shooting + Likelihood Maximization (LMax) approach [20, 21]. NC STATE UNIVERSITY Aimless Shooting and Likelihood Maximization • Aimless shooting generates an unbiased set of reactive trajectories by sampling the transition state ensemble. The procedure is as follows: (1) Find one point close to the transition state hypersurface, and build a reactive trajectory (i.e. one that connects the basins of reactants – A – and products – B – ) by shooting forward and backwards in time from that point (up to time T/2 and –T/2, where T is the total trajectory time). (2) Select, with equal probability, one of the points , , where ∆t is the timestep. Call this point the new initial (t = 0) point, and shoot forward and backwards in time again with a new set of velocities sampled from a Maxwell-Boltzmann distribution. (3) If the new trajectory connects A and B, accept it, otherwise reject it and start from the old trajectory. Go back to step 2. Store both the accepted and the rejected trajectories! 0x 0x 2t±∆x Aimless Shooting and Likelihood Maximization FIG. 2. Aimless shooting creates a sequence of interconnected trajectories. The three points on each trajectory are separated by a short time Af. = NC STATE UNIVERSITY 7 NC STATE UNIVERSITY Aimless Shooting and Likelihood Maximization • In summary, the algorithm to find a reaction coordinate using Aimless Shooting + LMax is: (1) Harvest independent trajectories using aimless shooting. (2) Propose a set of M relevant parameters, start using models for the reaction coordinate with one parameter (m = 1). (3) For each of the possible combinations of m parameters, maximize the log likelihood . Let be the maximum of all the values of obtained. (4) If of m = M, stop. If or m = 1, let m = m+1 and go back to step (3). • This process finds the best reaction coordinate for the given set of “important” parameters. RN m 1 1 ln 2m m R N+ − < 1 1 ln 2m m R N+ − > NC STATE UNIVERSITY Aimless Shooting and Likelihood Maximization • An example: polymorph transformation in terephtalic acid [22] NC STATE UNIVERSITY Outlook • Finding reaction mechanisms: – The Zero-Temperature String Method • Calculating reaction rates and defining reaction coordinates in complex systems: – Introduction – The Bennett-Chandler method – “Blue Moon” molecular dynamics – Transition path sampling – Aimless Shooting + Likelihood Maximization – Finite Temperature String Method – Other methods NC STATE UNIVERSITY Other methods • There are many other techniques beyond those discussed here to find rates or mechanisms of rare events, as this is currently a very active field of research. Some other methods are: – Parallel replica dynamics: This is useful mostly for “classical” MD simulations of complex systems. It consists on running a set of MD calculations in parallel until a rare event occurs. Then the simulation “clock” is advanced by the time taken by all simulations [15]. – Hyperdynamics: Similar in spirit to metadynamics and umbrella sampling, this method increases the energy of the reactants to reduce the activation barrier(s) [16] – Temperature-accelerated dynamics (TAD): As its name indicates, this involves running at higher temperature, where the rare event is not so rare. By assuming a form of the temperature dependence of the rate constant, this can be used to obtain the rate at lower temperatures [17]. NC STATE UNIVERSITY References [1] D. Masson, C. Richard, R. Martin, Int. J. Chem. Kinet. 8, 37 (1976). [2] D. Chandler, J. Chem. Phys. 68, 2959 (1978) [3] D. Chandler, Introduction to Modern Statistical Mechanics, Oxford University Press, Oxford (1987). [4] B.H. Mahan, J. Chem. Educ. 51, 709 (1974) [5] W.H. Miller, Acc. Chem. Res. 9, 306 (1976) [6] C.H. Turner, J.K. Brennan, J.K. Johnson, K.E. Gubbins, J. Chem. Phys. 116, 2138 (2002) [7] G. Ciccotti, M. Ferrario, J. Molec. Liquids 89, 1 (2000) [8] J.P. Ryckaert, G. Ciccotti, H.J.C. Berendsen, J. Comput. Phys. 23, 237 (1977) [9] H.C. Andersen, J. Comput. Phys. 52, 24 (1983) [10] See, for example, E.J. Meijer, M. Sprik, J. Am. Chem. Soc. 120, 6345 (1998) or B. Ensing, E.J. Meijer, B.E. Blöchl, E.J. Baerends, J. Phys. Chem. A 105, 3300 (2001) [11] M. Mugnai, G. Cardini, V. Schettino, J.Chem. Phys. 118, 2767 (2003) [12] C. Dellago, P.G. Bolhuis, P.L. Geissler, Adv. Chem. Phys. 123, 1 (2002) [13] D. Passerone, M. Parrinello, Phys. Rev. Lett. 87, #108302 (2001) [14] P.L. Geissler, C. Dellago, D. Chandler, J. Hutter, M. Parrinello, Science 291, 2121 (2001) [15] A.F. Voter, Phys. Rev. B 57, R13985 (1998) [16] A.F. Voter, J. Chem. Phys. 106, 4665 (1997); A.F. Voter, Phys. Rev. Lett. 78, 3908 (1997) [17] M.R. Sorensen, A.F. Voter, J. Chem. Phys. 112, 9599 (2000) [18] W.N. E, W.Q. Ren, E. Vanden-Eijnden, Phys. Rev. B, 66, 052301 (2002). [19] W.N. E, W.Q. Ren, E. Vanden-Eijnden, J. Chem. Phys. 126, 164103 (2007). [20] B. Peters and B.L. Trout, J. Chem. Phys. 125, 054108 (2006). [21] B. Peters, G.T. Beckham, B.L. Trout, J. Chem. Phys. 127, 034109 (2007). [22] G.T. Beckham, B. Peters, C. Starbuck, N. Variankaval, B.L. Trout, J. Am. Chem. Soc. 129, 4714 (2007). [23] W.E, W. Ren, J. Phys. Chem. B 109, 6688 (2005). [24] L. Maragliano, A. Fischer, E. Vanden-Eijnden, G. Ciccotti, J. Chem. Phys. 125, 024106 (2006)