 welcome to module 36 of chemical kinetics in transition state theory. This is a continuation of the last module. In the last module what we looked at is a calculation of transition state theory rate constant using molecular dynamics simulation. How we can apply molecular dynamics traits to do a thermodynamic integration and calculate rate constant? Today we will extend that a little bit and use the powers of MD to actually improve upon some of the assumptions of transition state theory as well. So, I will go back to a very old formula we wrote in module number 27. We can go back to that module and have a look. What we had really showed that the rate constant under the assumptions of a classical nuclei and Boltzmann distribution, okay thermal statistics being true can be written as this big integral. This integral I will just simplify just short hand notation for the purpose of this module. I will write this integral from dq2 to dq3n as dq over dividing surface. This notation we have used before as well in the course just to simplify our notation. This integral I will write as dp and there is no restriction here. It is over all momentum over all possible values of momentum. This I will write as e to the power of minus beta h transition state. Again that is a notation we have used before. H transition state really mean this h q1 dagger comma q2 to q3n comma p1 to p3n. So, this h is what I call h ts. Again the whole purpose is just so that I can write faster. If I write all of this all the time then we would not get anywhere. Multiplied by chi again I am dropping function of q comma p but that is always present. Just short hand notation and note that I am dividing by n. n is given here and so this is nothing but qr. Instead of qr let me write this integral clearly because this might get confusing since there is no h there. Integral dq integral dp e to the power of minus beta hr where hr represents the Hamiltonian of only the reactants as we have been using it so far. So, we use this formula in the last module to derive a slightly different formula for k ts t which is amenable to molecular dynamics. So, what we had shown is that we can play around with this formula little bit and write k ts t as this formula root kt over 2 pi m which is a thermal velocity into this function p and p is nothing but the probability distribution of being at q1. So, what we had discussed in the last module was that p of q1 naught is the probability density with q1 equal to q1 naught. So, if I set this one coordinate q1 as q1 naught and do a thermal integration over all other variables q2 to q3n p1 to p3n what will be the probability I will get that is what we defined as p of q1 and so this is p at q1 dagger where q1 dagger represents the transition state value and then we had worked out a scheme to calculate this p q1 dagger. We figured out that p of q1 dagger can be written as p of q1 r into this Boltzmann distribution where w is a kind of a free energy and then we figured out a scheme to calculate p of q1 r and this difference as well from last module. So, today let us do something more. So, this is our actual formula for k this big integral. Now, k TST is nothing but integral over this dividing surface dq integral over dp rho equilibrium into chi is replaced by this. So, we write this as a different form which is h of p1 I will define this in just a moment h of p1 equal to 1 if p1 greater than 0 this is equal to 0 if p1 less than 0. So, really it is a mathematical way of writing this line. So, if I write this h of p1 this integral will automatically be 0 if p1 is less than 0 and my integral over p1 by default will become integral from 0 to infinite which is what we have been doing. It is just an easier way of writing it I am just being a bit lazy that is all yeah. So, instead of a writing dp1 separately with a 0 to infinite I have written it in this simple functional form. If you want you can go ahead and write it in that particular form of 0 to infinity dp1 and dp2 to dpn what we have done in the last modules exactly the same thing all right. So, my point is this kappa this assumption that we have made today we are going to lift this assumption and we are going to use molecular dynamics to calculate this kappa more accurately than what transition state theory does. So, let us start I will simplify my life by defining something like this where kappa is nothing, but it is a definition. So, something stupid. So, if I divide these 2 dq dp rho equilibrium I will just write as e to the power of minus beta transition state chi p1 over m divided by dividing surface dq dp e to the power of minus beta h transition state h of p1 into p1 over m this is fine. But let us think of how we can calculate this expression have we made any real improvement we can keep on defining variables, but what about it how do I calculate this variable. So, let us give a little bit of moment in thinking what this chi really means chi is the probability of reaction if I am at that given q comma p. So, let us think about this coefficient a little bit more this chi with q1 equal to q1 dagger this chi tells me probability of reaction q equal to let me just say q naught and p naught equal to q naught and p equal to p naught. So, if my p and q are at this value here I start with this value on the transition state remember that transition state is simply specified by q1 equal to q1 dagger, but you have all these other variables with you as well which are q2 to q3 q3n and p1 to p3n they can be taking all kinds of values and depending what their value is your chi might not be 1. So, chi is a complex dynamical variable which depends on q2 to q3n comma p1 to p3n it is a function of it you give me these values and I should be able to figure out what this kappa is for you. So, how do we think about this effectively is accounting for recrossings. So, let us think about this statement a little bit what I mean by that if p1 is greater than 0 and I say chi is equal to 1 that is my transition state estimate. So, let us just come back to it. So, transition state theory makes a very simplifying assumption it says look only for the value of p1 forget all about your q2 to q3n and p1 to p3n p2 to p3n just look at p1 positive 1 negative 0. What is the underlying idea behind this assumption? So, let us say I have this potential barrier this is q1 it says that if I am moving in this forward direction I am reactive that is if I would have went back I would have ended in reactant if I move forward I will end in product. So, the reality is that we have a multi dimensional energy surface now. So, our energy surface let us say let me just draw 2 of them the kind of energy surfaces that we drew in the last module. So, let us draw something symmetric something simple this is let us say my transition state now a trajectory let us say comes from here this is let us say my reactant well it has a complex energy surface it sees it might turn around like this and come back. So, now when I am looking at let us say this point q1 comma q2 with the momentum that this point that the trajectory reaches this point with. So, this is my dividing surface this is the surface I am integrating over where I am setting the hitting setting the value of q12 now for this particular value of p1 comma p2 comma q2 my chi is basically 0 it is not reactive I can have another case as well. So, let us say I start here I come here I come here and I come here this is a bit more confusing. So, I have 3 points here that I am hitting 2 with positive momentum and it is reactive, but the problem is that I am counting this one trajectory 3 times when I do transition state theory. So, the main idea that we have is the following we want to calculate this kappa and I have written the full integral ones for you we want to do this full calculation we want to solve this big integral but using md. So, our idea is do constrained md with initial conditions q1 set at q1 dagger q2 to q3n p2 to p3n using Boltzmann distribution q3n p2 to p3n is equal to e to the power of minus beta h transition state. So, I sample all these q2 to q3n and p2 to p3n sorry p1 to p3n I am so sorry using Boltzmann distribution because that is what this factor I am getting here while q1 is set at q1 dagger and now basically we have to give a recipe of how to calculate this chi. As it turns out calculating chi using statistical mechanics is not easy there are attempts at it, but molecular dynamics is much more amenable to calculate this chi because this chi is truly a dynamical quantity or thermodynamics cannot really describe chi very well. So, let us look at a few different kind of trajectories a few different kind of crossings. So, I will making simple figures let us say on the x axis I have q1 and on the y axis let me draw time let us say this is my transition state geometry let us say at t equal to 0 my q1 is somewhere here which is a reactant and this side is product. So, q1 is my reaction coordinate the middle decides a transition state left of transition state is let us say my reactant right of transition state is my product good. So, as I move in time let us say the trajectory does something very simple. So, let us start with this one ok. This is the best case scenario for transition state theory it crossed the transition state exactly once. So, think about your transition state theory I am really sampling this big integral with q1 setting set to q1 dagger. So, this point here will also be sampled with q1 equal to q1 dagger and I although I have not drawn q2 to q3n and p1 to p3n they are having some value here and they have such values that I get this kind of a trajectory. If that happens my life is really really good my chi is simply one for this trajectory. Let us look at a more confusing trajectory. If this was the whole case then I do not need this module then kappa can be kappa is simply one transition state theory will be perfect, but that is not always the case unfortunately life is complex. So, there are trajectories the same let us say I start here let us say it hits here the same way, but here my values of p2 to p3n and q2 to q3n are different compared to this one. So, the values of other coordinates other than q1 at this point and this point are different that will lead to different dynamics now. So, this one let us say does the following just for example. So, now when I sample my e to the power of minus beta h transition state this point will appear this point will appear and this point will appear. I again emphasize this point corresponds to a certain value q1 is fixed, but it corresponds to certain values of q2 to q3n comma p2 to p3n. I will call this configuration configuration 1 I will call this configuration configuration 2 and I will call this configuration configuration 3. Now, note when I sample e to the power of minus beta hts the only constraint is q1 should be q1 dagger which is true for 1 2 and 3. So, all 3 points these 3 points will emerge in this distribution I will sample them with some probability, but they will be present. If that is the case and if I assign chi equal to 1 for this 1 let us say I use transition state theory basically. So, if I follow transition state theory chi equal to 1 for 1 equal to 0 for 2 equal to 1 for 3. So, transition state theory will essentially give me a total chi of 2. Let us assume that just for a moment that h transition state is has a potential has a same value at the crossings 1 2 and 3 that is not true, but for our discussion let us assume that. So, if I have that I got a total value of chi of 2 because transition state theory only looked at p1. So, p1 is positive here I am moving in the forward direction here p1 is negative here and again p1 is positive here. So, for this 1 trajectory I get chi equal to 2, but you see that is the problem that is precisely the problem that is the double counting problem of transition state theory or the recrossing problem of transition state theory. So, for this kind of a trajectory for this 1 trajectory I am getting chi equal to 2, but the correct answer is equal to 1 chi is really counting reactive trajectories that is a true meaning of chi. Chi should include all these three points as one point only because they are connected through this phase space they are part of one trajectory only. So, which means that transition state theory over estimates the rate constant by a factor of 2 in this case. So, a simple idea on how will I correct correction using md is start with p1 greater than 0 back evolve using molecular dynamics back evolve really means I am setting p to minus t remember that Newton's laws are perfectly time reversible. So, if you give me a value of q comma p I can also go back in time and tell you what q and p would have been sometime ago. So, that is a very easy to do in molecular dynamics forward evolve using molecular dynamics find number recrossings and set kappa equal to 1 over number of recrossings. So, let us just go by it in a moment ok. So, let us see what this will do has it solved anything I have just given you an algorithm, but does it solve this algorithm solve the problem or not. So, my sampling will only include this point and this point because I have said p1 is greater than 0 this one p1 is less than 0. So, this is not part of my sampling at all for these two points I will go back in time. So, I will let us say I start with 1 I go back in time I go forward in time I am sorry with p1 greater than 0. So, if I am setting at this point I will go backward in time I will go forward in time and I will do this trajectory basically and I will find that I have two points of recrossings one is this point itself, but by doing md I will also that there was a point 3. So, when I evolve forward in time from point 1 I will end up with point 3 which is also crossing in the with p1 greater than 0. So, for 1 my kappa will be set at half because toward crossing sorry forward crossing here gave me gives one recrossing. So, the starting point itself is a crossing plus when I move forward in time I found another crossing the total number of crossings became 2 for 3 similarly kappa is also half and it is half here because in backward direction now. So, if I start with 3 here now I will go back in time then I will find 1. So, when I propagate backwards in time I will have another recrossing. So, I basically I have 2 recrossings. So, now I have solved my problem. So, I get a half from 1 and half from 3 and that basically accounts for my factor of 2 here is accounted for. That is the main idea let me just add one more thing here. Let us imagine another kind of trajectory this is my time this is q 1 this is q 1 equal to q 1 dagger. Let us say start with a trajectory here and this trajectory let us say it does this this is not reactive at all for this trajectory chi must be equal to 0. I can have another trajectory let us say it does this for this as will chi should be 0. So, I this point will come this has p 1 greater than 0. So, I will be sampling this point I will be sampling this point. So, transition state theory would have said chi equal to 1 for this. Similarly, at this point transition state theory will put chi equal to 1, but this is wrong this is the correct answer. So, we have to account for this. So, the simplest idea is backwards evolution lead to reactants forward evolution otherwise. So, for this trajectory in my forward evolution I will end in reactant and hence chi will be said to 0 for this trajectory in the backward evolution I will end in products and hence chi will be 0. So, in short to calculate chi what we are doing a set q 1 equal to q 1 dagger and sample rest using e to the power of minus beta h distribution backward evolve forward evolve back evolution leads chi equal to simply 0. If forward evolution leads to reactants chi is 0 otherwise chi is equal to 1 over number of recrossings number of recrossings. So, this is the value of chi I will feed here and calculate these integrals. So, in summary we have today looked at how to calculate this transmission coefficient kappa using molecular dynamics to account for recrossings and this kappa will always be between 0 and 1. Thank you very much.