Sunday, October 25, 2009

Get a reaction: Intrinsic Reaction Coordinate


Lucas requested a post on computing an intrinsic reaction coordinate (IRC). An IRC is a special case of the minimum energy path (MEP) connecting reactant, transition state, and product (more about this at the end of the post).

There are many uses for the IRC, the most basic of which is to verify that the transition state you found actually connects the reactants and products you think it does (this is not always obvious from the normal mode associated with the imaginary frequency). It is also a good way to identify any reaction intermediates you haven't thought about. The IRC is also a way generate an animation of the reaction (see for example ChemTube3D).

An IRC is generated by following the gradient down-hill from the transition state to products and reactants. (It's a bit like a geometry optimization initiated at the transition state, but the algorithm used to follow the gradient is different). Computing a complete IRC thus requires a minimum of two separate runs, and you must find the TS before you can compute an IRC.

The screencast shows how to generate the IRC associated with the 3TSa transition state I found in a previous post. The first input file I generated contain the 3TSa coordinates and the following keywords:

$contrl runtyp=irc $end
$basis gbasis=pm3 $end
$irc saddle=.t. tsengy=.t. forwrd=.t. npoint=100 stride=0.2 opttol=0.0005 $end


saddle=.t. tells GAMESS that the coordinates in $data corresponds to a TS. Because the TS is a stationary point the gradient is zero (or very small), so there is no gradient to follow down-hill. Therefore, the IRC is initiated by displacing the geometry along the normal mode corresponding to the imaginary frequency. Therefore, saddle=.t. GAMESS expects a $HESS group, which contains the frequency information.

Whether the normal mode points towards the reactants or products is completely arbitrary, and the forwrd=.f. keyword is used to change the direction in which the normal mode points. forwrd=.t. is the default and I only include it in the first run for completeness.

tsengy=.t. tells GAMESS to re-compute the energy of the TS (re-compute because you already computed it when you computed the frequencies). If you want to make an energy plot in MacMolPlt as I show in the screen cast, you must set tsengy=.t.

The displaced geometry will have a non-zero gradient, the negative of which points downhill. The geometry is displaced in this direction, the gradient is recomputed, and the procedure is repeated npoint=100 times or until the (root-mean-square) of the gradient is below opttol=0.0005 atomic units, at which point you are at a minimum. (Notice that the default for npoint is 1, and must always be changed).

stride=0.2 tells GAMESS how far along the gradient it should displace at each step (in this case 0.2 atomic units). The default is 0.3, but when I tried this, the IRC calculation in the forward direction failed after 7 steps, because the coordinate deviated too much from the minimum energy path, so I had to decrease it. This means I will have to take more steps to reach the minimum.

As you saw in the screencast, GAMESS takes 100 steps in each direction, which means the IRC does not go all the way to the minimum on each side of the TS. GAMESS does print out restart information to continue the IRC, which I didn't do in this case.

Finally a note on MacMolPlt. As I show in the screencast, MacMolPlt can be used to combine the two IRC runs (forwrd=.t. and forwrd=.f.) by using Files > Add Frames from File... You should first load the file that goes to product and then add the file from the IRC that goes to reactants. Here you should tell MacMolPlt to reverse the order of the structures in the file, so the animation goes from reactions to transition state to product.

This is done by selecting "make these points negative". "Negative" refers to the x-coordinate in the energy plot. The x-coordinate is what separates an IRC from a minimum energy path. An IRC is a plot of the energy of each structure along the minimum energy path vs the root-mean-square change in the mass weighted Cartesian coordinates relative to the TS structure. The TS thus has a x-coordinate of zero, while structures leading to the reactants are defined by negative RMSD values.

If the energy is plotted against something else, such as a bond-length, then it is simply a minimum energy path. Because an IRC uses mass-weighted coordinates, using different isotopes will lead to different IRCs, and one use of IRC is to determine tunneling corrections to the reaction rate.

27 comments:

Unknown said...

thanks jan,
and sorry for any disturbance.

NUchem said...

Jan,

Thanks also. I've always wanted to know how an IRC was used. Great job.

NUchem

Jan Jensen said...

Lucas - no problem. Thanks for the suggestion

NUchem - Thanks! Glad you liked it.

Alonso Rosas said...

Hi!

I am a student of Chemistry and I have found very useful your blog. I have just a question about some of your post. I am very interested in jmol and I have read about but I have not found yet how to install jmol in a mac computer. I hope you can help me.

Thanks!

ALonso

Jan Jensen said...

Alonso - downloading and installing jmol is described here here. Hope that helps. If not let me know.

Alonso Rosas said...

Installing jmol can't be easier!

Thanks and sorry for any disturbance

ALonso

shuai said...

hey jan
your blog is so amazing
i am a student who is learning how to use gamess and your blog really answer me lots of questions
i have a question
if I want to optimize more than two water moleculars (like three water molecular in constraints such as the distance between each O is 2A), how to set up the NONVDW and IFZAMT?
thankyou
shuai sun

Jan Jensen said...

Alonso - glad to hear it worked

Shuai - I am very happy to hear you like the blog.

Now to your question. Let's say the oxygens are numbered 1, 4, and 7:

$zmat nonvdw(1)=1,4, 1,7, 4,7 ifzmat(1)=1,1,4, 1,1,7, 1,4,7 fvalue(1)=2.0, 2.0, 2.0 $end

shuai said...

hey Jan
thankyou for your reply
i tried three water molecular optimization and it worked!!! thankyou
then i tried to optimize a water molecular surrounding by two CH3CH2COO- , and i replaced two of the water oxygen atoms in ifzmat to the oxygen atoms of O-in acid molecular, but gamess told me that SCF cannot be converged
do you know why?
thankyou
shuai sun

Jan Jensen said...

Shuai - I suggest posting this question to the google gamess group.

When you do, be sure to include a copy of your input file.

Shuai said...

thankyou
people there helped me to figure it out
thankyou very much
Shuai

shialchemy said...

Hi Jan,
Usually I carry out IRC calculations to verify the saddle point is the transition state of interest.
Have you ever met the situation that the energy of structure along the reverse MEP is "lower" than that of reactants?
I am dealing with bi-molecular reactions. I guess it is likely that the complex formed by two reactants have lower energy than the separate reactants. Interestingly, in my case, B3LYP gives reasonable IRC plot (no negative energy compared to reactants) but M06-2X gives complex with lower energy.
Do you have any comments?

Jan Jensen said...

Molecules always attract in the gas phase, meaning they always form a complex that is lower in energy than the separated molecules.

I have written about a related topic here: http://molecularmodelingbasics.blogspot.com/2012/01/negative-activations-energies-in.html

Jan Jensen said...

PS: just to be clear, I am talking about the PES here, not the free energy surface.

shialchemy said...

Thanks a lot. I have seen your suggested post. Very helpful.
I had another question following that post. You suggested that if we are modeling a reaction in solution, then the energy barrier is E(Ts)-E(R1/R2).However, in most cases, we are modeling reaction in gas phase but compare the rate constant and energy barrier measured in solution. Then which one is a more reliable comparison, E(Ts)-E(R1)-E(R2) or E(Ts)-E(R1/R2)?
I know that I am not performing an apple to apple comparison but solvent effect on energy barrier is so difficult to study.

Jan Jensen said...

E(Ts)-E(R1/R2)

The usual caveat applies: gas phase calculations on very polar molecules in highly polar solvent are poor approximations.

shialchemy said...

Sorry to bother again. This discussion is very beneficial for me, but I still have some issues about the correct way to calculate the barrier.

It seems to me that E(Ts)-E(R1/R2) is a more physical description for bi-molecular reaction, as two reactants should form some kind of complex before reaction, no matter in gas or solution.

But if we use E(Ts)-E(R1/R2) to calculate the energy barrier, does that mean the reaction now becomes a single molecule (first-order) reaction? Of couse, I am talking about the rate constant obtained using transition-state theory.

Does this suggest there is no "real" bi-molecular reaction?

Ideally, such problem can be resolved if we can perform first-principles calculation for reactants and its environment together. However, as we both know, it is impossible now, not even in near future. And it seems that people usually just optimize the reactants separately and somehow ignore the complexion energy.

I am not sure whether which one is a physically right way, especially when in current stage, we usually perform gas-phase calculation

Jan Jensen said...
This comment has been removed by the author.
shialchemy said...

I hope my long comment did not annoy you. Now I think entropy also plays a role. Generally speaking, if there is no strong intermolecular interaction, such as hydrogen bond, it is always to safe to optimize reactants separately without considering the formation of complex in gas phase calculation.

Jan Jensen said...

Not at all. I have been thinking about the problem, but was too busy to do too much else.

The bottom line is that I have changed by mind, and updated the post accordingly: the activation energy for a bimolecular reaction should generally be computed as E(TS)-E(R1)-E(R2) - also for condensed phase reactions.

There are probably exception to this in the condensed phase, but this is the best rule of thumb I can come up with.

rfpeters said...

Hi,

Absolutely wonderful blog! I have been using GAMESS for calculations of Raman intensities. Have you done any work with Raman calculations? My results don't match too well with my experimental results. I thought I would check to see if you have any tips.

Jan Jensen said...

I'm very happy to hear you like the blog!

I haven't done any work on Raman intensities, so I can't really help you there. I suggest searching the literature or posting your question to ccl.net

Anonymous said...

Hi Dr Jensen

Just got my first IRC to an extremely large part through following your blog.

You describe combining the two files (forward and back) using MacMolPlt ... how is this process done in the absence of MacMolPlt?

I put my structures on an intranet and rely on jmol.

In anticipation, thanks ...

Jan Jensen said...

Right off-hand I don't know of another code that does this. I am not sure, but perhaps you could still combine the files with MacMolPlt and then save the file for use with Jmol. A

Alternatively, you'll need to write a little script that reverses the order of the structures for the backwards IRC.

Anonymous said...

I used MacMolPlt - it couldn't have been easier, thanks.

Anonymous said...

Thank you soooooooooooooooooooo much.
Your blog is very very very useful for me.

Jan Jensen said...

Thank you! Very happe to hear it.