Hey, everyone. It's Corin here from Rowan. In this video, we're going to look at how we can go from a simple 2D drawing of a given reaction to a prediction of what the barrier height for this process will be. In other words, what is the energy difference between the starting materials and the transition state, and what will the rate of the reaction be? This is a pretty useful computational task. It's often one of the ways computation is introduced, or ways students might see computational values being used.
But it's actually somewhat complex to do and it requires a couple of different steps. So in this video, we're just going to walk through it from start to finish. And by the end of this video, we'll be able to predict what the barrier height between cyclopentadiene and maleic anhydride, how fast this Diels–Alder reaction will happen. So this is just sort of an illustration on Master Organic Chemistry. And we're going to take this.
Let's pop over to Rowan and start drawing these materials in. So we'll go to our tutorials project and we'll make a new folder for this video, which we'll call maleic-nehidride plus cyclopentadiene. We'll create this and we'll go in here so all our calculations are in one place. Okay, so to model this process, we're first going to need to model the starting materials and then we're going to need to model the transition state.
So let's start with the starting materials, which will be a little bit simpler. We can navigate to the calculation workflow and we can draw the man in 2D. So our first starting material is cyclopentadiene, this five carbon cyclic diene. What we can do is go down here. There's actually a template just for this exact species. So we'll say save. We'll convert it to 3D.
And for this video, let's do everything at the UMA small level of theory. So this is the "universal model of atoms" from Meta's FAIR-chem team and a bunch of other folks. This model was published earlier this year and is a good all around choice for a lot of chemical tasks. We'll leave this in the gas phase, although we could select a solvent if we wanted to, but we'll just do it in the gas phase for now. And then we'll select "optimize" and we'll select "frequencies."
So we'll say, let's rename this so it's a little more human readable. We'll say cyclopentadiene and we'll say "submit calculation." This calculation should start running pretty quickly. So now our other reactive partner is maleic anhydride. So this will be our dienophile. And so we can do the same thing here. We'll go to the 2D structure editor. We'll select this template. We'll type in our molecule.
We can just double check we've got the structure right. We do. Now we'll convert this to 3D. We'll call this maleic anhydride. We'll say UMA here—very important to double check that we run every calculation for this video at the same level of theory. And then we'll run this for optimized frequencies as well.
OK, so let's check back on our pentadiene molecule here. And yeah, mean, frankly, this looks pretty good. This looks like it's progressing OK. So now we have to actually find the transition state for these two species together. And while we could do this a number of ways, so we could just start a transition state optimization here just from a guest transition state, or we can conduct a scan, probably the best way to do this is to use Rowan's double-ended transition state search functionality, which we released just a few weeks ago.
And then at this what we'll do is we'll have a reactant and we'll have a product and we will sort of like automatically find the transition state between them using the frozen string method. So to get started, let's just draw the product we want, which is this molecule here, because this will actually be a little bit easier, I think, to draw than the reactant complex. It'll be easier to go from the product to the reactants than from the reactants to the product.
Really, we should be drawing this in 2D. So let's start there. If we look here, we have a six-membered ring with a single bridging methylene here. And then we have the five-membered ring. So just as before, whoops. We can do this. Okay, I want to drag this around here. Take this and align those two, then make this a double bond. And let's check that this matches. So it does. We'll save this and then we'll convert it to 3D.
And so this now we can see is, you know, if we scroll down a little bit, this half corresponds to the maleic anhydride portion and this half corresponds to the cyclopentadiene portion. So we'll save this as our product. Then we'll copy the XYZ coordinates and it'll go up here. And we will say "paste XYZ" here. So we'll submit. And now what we need to do is we need to actually separate these apart into a reactive complex. And there's a couple of ways we can think about doing this, but perhaps the default is just to remove the bonds. So we'll click on this bond, go to default and we'll say remove bond. Then we'll just double click on this structure, hit Z and drag it away. And there we go, two separate species.
And one thing we can do even just to clean this up beforehand is say optimize with GFNFF. That will just sort of let these molecules relax a little bit so that the optimization is simpler later on. Okay, and so we'll call this "double ended TS search." We'll say optimize inputs, optimize transition state. We'll select UMA and then we'll leave this in the gas phase. And what this is going to do is it's going to try to find a transition state between our two starting materials and our product. So one really important thing to note about this is that it takes one single structure that has both species in it. It doesn't take two separate structures. So we have to have actually drawn them in the same canvas here. So we'll say submit and this will start running shortly.
Okay, so this is our converged structure for cyclopentadiene. We can see that all of the frequencies are positive, which means that this is actually a minimum and not a saddle point. And we can say the same thing is true for maleic anhydride. All the frequencies are positive. Now the key tab for the purpose of this video is actually this thermochemistry tab. So here we have the electronic energy, which is sort of the default energy that you see, for instance, in an optimization. We also, since we ran a frequency calculation, have the total enthalpy and the total Gibbs free energy. And this we get by adding these sort of vibrational, translational, and rotational corrections to the molecule. So this more closely represents sort of the free energy ah with some entropic components. So we have some entropy in the vibrational modes. We have some rotational entropy. We have some translational entropy. There's obviously approximations behind all of this. But if you want to actually predict how fast a reaction will happen, you want to be using Gibbs free energy.
And so let's go over here and let's actually start making a table with all of these values pulled in from these structures. We'll say, okay, we'll say a cyclopentadiene, say maleic anhydride, then we'll put room for the transition state here. And let's just zoom in a little bit for the sake of this video.
All right. One thing I actually should do is insert a column to the left so we can say what this is. So do it. Energy, enthalpy and Gibbs free energy. We'll start pasting values in. So in Rome, you can just click and it'll automatically copy it to the clipboard. So for cyclopentadiene, the energy is -194. The enthalpy is also about -194. And the Gibbs free energy is also about that. And so these values right now are in Hartrees.
So a Hartree is a very, very large unit when you're thinking about differences in energy between similar structures. So for a difference between two conformers or between a ground state and a transition state, a Hartree is a lot of energy, right? So one Hartree is 627 kcals per mole. And the reason these numbers are so big is that this is actually the interaction energy between every nucleus and every electron, relative to nuclei and electrons at infinite separation. And so we expect these values to be really, really large. And if they're not, that usually means something's wrong.
All right, we can do the same thing here for maleic anhydride. So we'll take electronic energy, put it here, enthalpy here, and Gibbs free energy here. All right, we're just going to wait one sec for this double ended transition state search to run and we'll fast forward to witness running.
All right, so I fast forward just a little bit to where the double ended transition state search workflow has started running, just so we don't have to stare at a blank screen for a few minutes. And we can see that what we're doing is we're starting over here. So these structures are reactant structures, and then these structures are product structures, if we go all the way over here.
And now what the workflow is doing is it's finding the string that connects these two structures. And so ideally, what we'll see through here is that we'll see that energy goes up. And then it starts to go back down and sort of at the point where it is up will be a good guess for the transition state. And this is sort of just the classic way the double ended transition state search method works.
So one other thing that we could do, I just realized, is that if we really care about what the overall energy change in the reaction is, we can actually resubmit this structure and look at it as the product. So let's do that now just for illustrative purposes. So let's say here, optimized frequencies. We'll call this just "product," because I don't want to try to systematically name this on the fly. And then we'll see again that this is using the UMA small model, which is great.
OK, oh, here we go. So now we're finding something that starts to look a lot more like a transition state. Yeah, this is looking quite good. So we can see here, yeah, starting to bend over. And there's still this, you know, the workflow is not quite connected, but these structures here will definitely be pretty good guesses for their transition state. And it's likely that this will converge when it automatically re-optimizes the transition state next. We can also always go back and check on the product just to make sure that this is running okay. And I bet we'll start to see some optimization steps coming through shortly.
Alright, so we see that this just finished. And now not only is the entire string complete so we can run through and watch the reaction happen, it's also automatically optimized a guest transition state so we can look at this. Sort of see what this is. We can view the underlying calculation. Yeah, we can look at the frequencies and say yes, there's sort of a single imaginary frequency here.
This looks quite good. And then we can see that there's thermochemistry parameters here as well. And so this is for the transition state.So let's take this and start pasting this into our Diels–Alder table. So we'll take the electronic energy here. We'll take the enthalpy here. And we'll take the Gibbs free energy here. And then one last thing, because I think I saw that it finished, is the product. And so here again, frequencies are all positive, and we'll take the thermochemistry values from the product, and we'll put them here. So, just a little tedious here. Sorry about that. Maybe in the future we'll have a better way to do this.
Okay, so we finally have all of the energies we want. And now the question is, how do we turn these numbers, which as you recall are all sort of like ludicrously large energy values that correspond just to infinitely separated nuclei and electrons, and how do we convert this into a potential energy surface? And what we basically need to do is we need to sum up the energy of each of these and then look at the difference in energy.
So there's a few ways we can do this. In this case, I'm just going to sort of say, cyclopentadiene, so Cp plus Ma, maleic anhydride. And we're just going to add these two energies together. And what this does is this sort of says, what's the joint energy of the system if it's not interacting at all? And what we can see now is sort of that we have the reactants together. We have the transition state and we have the product. And so these now, you can tell these numbers are all quite similar. And that's a good sign. That means that the differences in energy between them are going to be pretty small, like we expect.
So we can sort of create a new thing, which is delta E, delta H, and delta G down here. And now we'll look at the relative energies. So we'll say, these are this minus this. And we'll say, since this is in arteries, we'll want to convert to kcal per mole. And we'll do this by 627.509. Okay, and now we'll want to do the same thing over here as well. So if I were slightly better with spreadsheet formulas, I could just copy it automatically, but I didn't freeze the column.
Okay, so now what we have is energies in kcal per mole that are relative to sort of the zero point here. So this would be zero in each of these cases. And so we can say that in electronic energy terms, the transition state is 10.8 kcal above sort of the cyclopentadiene plus maleic anhydride separated. And then the product is 33 kcal per mole downhill. But this is missing something important. Remember that this doesn't take into account entropy at all.
And so if we sort of take it like this and we just try to predict the barrier this way, we're getting an electronic energy barrier. And this doesn't, this is just wrong, right? There's an entropic cost to bringing the two molecules together into a complex of sort of depriving them of that translational and rotational freedom. And that's only properly taken into account if we look at the Gibbs free energy, which is this Delta G. And so here we can see that the Gibbs free energy barrier, the correct barrier for reactions that we're observing under experimental conditions is actually 23.8 kcals per mole. And so this is now, instead of a reaction that would be super, super fast at room temperature, this reaction would probably require a little heating, you know, if it truly is a 23.8 kcals per mole barrier. And now the reaction is still downhill, but it's not quite as downhill because there's sort of an entropic cost to going from two reactant molecules to one product molecule.
Okay, so if we were reporting something in a paper, we would say that UMA predicts the barrier to this reaction to be 23.8 kcals per mole. And we can sort of, you know, fix the sig figs here so it's less uh ludicrous. Just one more point I want to illustrate here, which is sort of, think, something that people often get wrong. That's really quite bad. And that's if we take this system here and we resubmit it. um We'll say we'll just resubmit for a current structure. We'll go here. We'll say pre-reactive complex, say, "optimize" + "frequencies," UMA, that this is not a good structure to calculate a barrier from. We need to calculate the barrier from two separate species and then add up the free energies of them separately, because barriers from pre-reactive complexes are almost always significantly lower than the barriers from isolated reactants.
And that's because when we form this pre-reactive complex, we haven't taken into account the entropic cost of bringing these two species together. And so sort of once they're already in close proximity to each other, the reaction might be rather fast. But you know, back of the envelope, it takes about 10 kcals per mole to bring two species together. And that's sort of priced in here, even if they're formally non-interacting in an electronic energy sense.
Okay, and so if we actually look at this optimization, we can sort of look at the frequencies here and look at the thermochemistry. We go and take these values—insert one column left so as a CP plus MA one comp—what we can do is we can take the electronic energy and we see that it's very similar. We can look at the enthalpy. We can look at the Gibbs free energy here. And what we can actually do is do the same operation so we can do to say this minus this times 627.509.
What we see is that although this complex is allegedly, you know, 4 kcal per mole downhill on the electronic energy surface, it's actually 7 kcal per mole above the separated species when we take into account free energies. And so if we were actually looking at the reaction barrier relative to this complex, which we see isn't actually, you know, properly uh a true ground state complex. So we need to, you know, it's not clear what exactly the minima of this pre-reactive complex is, but just sort of setting that aside for now, we would actually be underestimating the barrier of the reaction by 7 kcals per mole. So this is sort of a common mistake, and I'm here to tell you not to do it.
Anyhow, hopefully this was a useful tutorial for how to find transition states and barrier heights. And I wish you all the best of luck in doing this for your own systems.