 Welcome to module 35 of chemical kinetics and transition state theory. Today's contents are going to be somewhat different than what you will typically study in your courses. But I believe it is still very important to learn this point. Remember we are learning so that we can calculate rate constants for a general problem. What we have learned so far is we started with these rate theories. We looked at transition state theory in great detail. And in the last few modules, we looked at molecular dynamics as an alternate way to calculate rate constants. In the last module, we looked at how rolling spheres on potential energy surfaces can lead to more insights on chemical reactions. Today we will ask a different question. First thing is if I can run these molecular dynamics. Are rate theories really necessary? Well, they are. The reason is these molecular dynamics are very computationally expensive. Specifically, if rate constants are slow, if you let us say have a reaction that is happening on a millisecond time scale, running a molecular dynamic simulation for that longer time is extremely hard. Very few computers in the world can do that kind of simulation. And it is not needed necessarily. We have transition state theory or alternate rate theories if you wish that can calculate these rate constants anyway. So, if all that you want is to calculate a rate constant, then MD is not necessary. It is one of the tools. So, transition state theory is computationally cheap. It is effective. However, and we looked at how to use transition state theory to actually calculate rate constants as well. We can divide partition functions into translational, rotation and vibrational and then calculate a rate constant. However, there are situations when this separation of partition functions is not always possible. It is not always accurate. So, is there a way for me to use the power of molecular dynamics simulations to calculate transition state rate? Actually answer is of course, yes. And it is commonly being used. So, I wanted to give you a little glimpse of it in this module and the next module on how we can use the power of simulations to calculate a transition state rate. So, we will have to go back to a rather old expression I derived some time ago, almost 6, 7 modules ago. We derived transition state theory in two ways. And in the second approach that we looked at, this integral we had derived. If you do not remember, you can go back to your module 27. This is nothing but this integral is your forward flux. This is an integral over the dividing surface of e to the power of minus beta h. And this is nothing but your partition function of the reactant. So, this is going to be my starting point. I will manipulate this thing a little bit more. So, let me start with this. So, I will do the following. First thing, I will separate out all momentum integrals completely. In the denominator, I will take out the momentum term. So, I have taken out this term and h r I am saying is equal to p i square over 2 n plus v of r multiplied by all kind of other integrals d p 2. I am also assuming this h Ts is equal to sum over i p i square over 2 n plus v of transition state q comma. So, I will have to delete this thing because I will run out of space. I will just write it below for your convenience. And in the denominator, I will have the actually the exact same term, both of them having the same limits. You will notice I will have this for all 3 n integrals and then I will have the position integrals that I cannot do anything about, not yet. I will simplify this a little bit. What does this mean really? This term is really the potential of transition state at q comma p really means v of the reaction coordinate equal to q 1 dagger the transition state geometry comma all other coordinates. This is basically the definition of v PS that my v is set of q 1 is set at q 1 dagger which is my transition state geometry. So, this is q 1 equal to q 1 dagger comma q 2 to q n q 3 n sorry divided by an integral over q 1 to q 3 n e to the power of minus beta v q 1 to q 3 n. So, I have just written everything very very clearly every term explicitly. So, all the momentum terms here cancel. Let me just simplify this term. So, this term that I have well we have been doing this kind of integrals very frequently. Today I have forgotten to provide you the integral equations here, but that is not a big deal. These you can look up your integration tables. The numerator is k B T and the denominator is going to be root 2 pi m k B T. So, these integrals are just looking up integral tables. Do not worry too much about them. So, this is nothing but root of k T over 2 pi m. So, that is nothing but your average velocity average momentum. So, I end up with this equation now. The all the momentum integrals had cancelled except for p 1 and the p 1 integral gave me this term k T over 2 pi m square root and I am left with still the integrals over all other coordinates d q 2 to d q 3 n d q 1 to d q 3 n this this. Now, the fun will begin. I am going to define a function p of q 1 not to be this where instead of q 1 dagger I have written q 1 not here. So, just be careful here I have this is not a mistake I have consciously written this. So, q 1 not is some variable. It can take any value not necessarily only the transition state value some q 1 not call it x if you want call it y if you want I will call it q 1 not because that is what I want and the denominator is exactly the same. So, k is nothing but this square root of k T over 2 pi m term into p of q 1 dagger specifically. So, if I put q 1 dagger here I will get q 1 dagger here and I will get exactly this term sound something trivial to do something stupid to do, but let us do it nonetheless. Let us now think of what I have defined here this p of q 1 not actually this p of q 1 not is nothing but average of q 1 minus q 1 not delta. So, let me just be clear average of any quantity o over reactant this I am defining to be 3 n divided by sorry operator o here e to the power of minus beta v r because that is how averages are defined that is how thermal averages are taken and here I am specifically taking thermal average with a denominator specified to only reactants that is how transition state theory does to you. So, I am saying that this p q 1 not is nothing but average of this function and I do not want to get into the technical nuances here this is called Dirac delta function, but if you do not know then no problem physically this is p of q 1 not is the probability thermal probability of finding q 1 equal to q 1 not. So, if I take this integral here and I think that q 1 is fixed there is no integral over q 1 q 1 is a specified q 1 is specified to q 1 not in the numerator then I will get exactly this this thing I do not worry about this I have included in my potential. So, do not worry about this part that much. So, p q 1 not effectively tells me if I am at thermal equilibrium what is the probability that q 1 is equal to q 1 not particularly it is a probability density to be more accurate. This p q 1 not is there any way we can calculate it. So, that is the question remember a rate constant my transition state theory rate constant is nothing but some constant that I anyway know I know temperature I know masses and all this into p. So, if I can tell you how to calculate p I can tell you how to calculate rate constant. So, actually molecular dynamics can calculate p that is the point. So, first thing is what was molecular dynamics again molecular dynamics was xi dot equal to vi vi dot equal to my acceleration which is nothing but minus 1 over mass well this thing the solving these equations gives you NVE ensemble. But actually there are ways of solving equations that gives me NVT ensemble as well. I am not going into these details, but effectively you can add a thermal bath to your simulation it is at the same way as Brownian dynamics work you think of this system which is open which has been constantly being colliding from other particles at some temperature and that can be simulated there are many many different ways of doing that. So, for now just assume that I can simulate an NV2 ensemble as well. Well pq1 not then it is not that hard to find really pq1 not I can find using md using the following way run md simulation this NVT ensemble. So, you have some energy surface with you I do not even care about this energy surfaces I am trying to ask you what is the probability density of finding the particle at some q1 not well the idea is I define a very small distance here I asked so choose a small delta q. So, in my md simulation the point is how many times ratio of times that q1 was equal to q1 not divided by total simulation time. So, what I am saying is if I run let us say a 1 nanosecond simulation and out of that 1 nanosecond q1 was equal to my system was exploring this region for 1 picosecond then the probability that I am at that region is 1 picosecond divided by 1 nanosecond right or I explore this whole region in 1 nanosecond. And in that 1 nanosecond I was here only for this much time only for 1 picosecond only if I take the ratio I will get the probability density at that point. Well there is one little trick here in an md simulation you have discrete times so q1 equal to q0 mathematically might not happen you might come very close to q1 not but not exactly on top and that is a little just nuance and the way to solve that nuance is that we choose this small dt delta q and so we find this ratio. So, time when q1 is really close to q1 not instead of equal to q1 not so when q1 minus q1 not is less than delta q. So, q1 minus q0 q1 not is very small divided by simulation time and since I have included this delta q I have to divide by delta q because if delta q is large then well you will spend more time there. So, if I instead of choosing this box if I choose this box well of course you are going to spend relatively more time so we divide by delta q. So, this is the prescription using md you can actually find p of q1 not now and once you have p of q1 not you can also calculate p of q1 dagger and remember the transition state rate is proportional to this term multiplied by some root of kT over 2 pi m but we actually have another snag here another problem here. The problem is if the barrier height is quite large if this is much larger than kT then if you run an md simulation the time that you will spend here is extremely small most of the times you will be exploring this region or you are going to explore this region and once in a blue moon you will actually explore this region. So, if I run a 1 nanosecond simulation you might spend 1 femtosecond there so it is not very efficient we have not achieved much here so far we have to just hit this remove this little trouble that we have and there is a very clever way of doing it so I just wanted to show that and this is a very common trick we define a function now just bear with me you see this is going to be very beautiful minus kT ln of p of q1 not just bear with me and let me define this function for those who are familiar with a little bit more thermodynamics and statistical mechanics this is a definition of free energy this is how free energies are defined in stat mech but do not worry if you have not seen this you can take this as a definition all right this will give me p of q1 not to be e to the power of minus w of q1 not divided by kT fine that will give me p of q1 dagger to be e to the power of minus w q1 dagger over kT fine let me also define p of q1 reactant so I am thinking of this energy surface this is q1 dagger and let us say this is q1 r something that is close to a minima fine seems like I am doing some very stupid things I am trying to confuse you or fool you but I promise you I am not or I will take a divide these two so p of q1 dagger equal to p of q1 r into e to the power of minus w of q1 dagger minus w of q1 r over kT I just taken a ratio of these two terms still seems silly I have not achieved anything looks like my point is I will now figure out how to calculate this and this separately p of q1 r is easy q1 r is close to a minima so I can actually use md effectively to calculate p of q1 r why because your simulation is going to spend a lot of time close to the minima so if I run a 1 nanosecond simulation for more than half a nanosecond you will be close to a minima so you can get good statistics you can quickly calculate p of q1 r because trajectory spend more time close to q1 r good but what about that I still have to calculate this w and w is some weird the minus kT ln of whatever so now let us see what we can do this is a real brilliance it was originally given by Eyring and Keck Poliani these people developed this I have to calculate this quantity and Bennett and Chandler eventually made it very concrete in terms of a simulation I am going to write this difference as an integral you are rarely going to see this but I take a simple subtraction of two terms and convert it into a complex integration but it helps okay so you can quickly verify that this equation is true if I take the derivative and take this integral I will simply get w w at q1 dagger minus w of q1 r fine now the point is this del w over del q1 not let us try to calculate that so w was defined earlier like this here and p I had already defined like this if I want to take del w over del q1 not this is nothing but minus kT into 1 over p of q1 not into del p over del q1 not so I take the look at this p q1 not is only here so this is equal to minus kT over p of q1 not and I write all these integrals now and I take the derivative of this term the derivative of this term is e to the power of minus beta v q1 not q2 q3 n into derivative of the exponential so just differentiating by parts nothing fancy minus beta del v over del q1 not divided by the denominator is independent of q1 not so I leave the denominator alone now notice beta is nothing but 1 over kT for a minus here and a minus here becomes plus I have kT into beta divided by p so I take this p and write it again so I get let me write this term first sorry for a little confusing statements q3n into del v over del q1 not divided by you notice once I divide by the p here if I am going to divide by p like this this denominator is going to cancel this integral dq1 to dq3n is exactly written here so this exactly cancels and I am left with and I notice kT into beta is of course 1 beta is 1 over kT ok kT into beta is 1 ok so I am left with this integral so this integral is actually called that is w is called the potential of mean force what you have gotten is this is nothing but average of del v over del q1 not at q1 equal to q1 not so I am taking a thermal average and I am fixing q1 equal to q1 not both in numerator and denominator all right del v over del q1 not is nothing but minus of force so dw over dq1 not is nothing but average thermal force at q1 equal to q1 not so that is why it is called potential of mean force you are taking a mean force and you are converting a potential out of it ok ok so this dw over dq1 not actually can also be computed using an MD simulation ok so the point is in MD simulation there is another trick I can do that I am not covering in this course I can put constraints constraint q1 equal to q1 not use NVT ensemble so if you download any of these popular software like gromax or lamps which can run MD they all will be able to do these tricks for you ok so once that happens you can basically run MD simulation you have constrained q1 equal to q not and that can be any value that can also be the transition state value and once you have done that you can calculate average force so you calculate this derivative with respect to q1 run your simulation and fix q1 not q1 equal to q1 not and just run the simulation and just find this average and that will give you dw over dq1 not ok so now we have a full prescription pq1 dagger we had shown was pq1r into this exponential and we had written this exponential as an integral so this I have written as pq1r into e to the power of minus beta integral from q1r to q1 dagger del w over del q1 not dq1 not ok pq1r I can calculate directly from an MD simulation by using this kind of a trick and I do this integration numerically I start at q1 q1 not equal to q1r I calculate this value I shift q1 to q1r plus some dq I recalculate this so this integral itself now can be calculated numerically at using small grids of q1 so what I am again doing is I am setting q1 equal to q1r here I calculate dw over del q0 then I shift q equal to some q1r plus some del q and I recalculate this function here and I keep on calculating it till q1 dagger and this integration then is nothing but just f into del q sum over f into del q so this integral is nothing but sum over i fi del q and again this dw over del dq I can calculate using a constrained or NVT MD simulation so this was a quick summary today of how we can use an MD simulation to calculate a transition state theory rate so this is useful when partition functions are not easily available when you cannot separate let us say rotation and vibration or if the vibrations are not is harmonic then what do you do when you do this kind of potential of give force thank you very much