Sunday, October 25, 2009

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:

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.

The movie (which is best viewed on the vimeo site using full screen) very elegantly shows the evolution of the program (notice the time on the lower right hand side of the screen), from the very beginnings in 2006 when Geoff and Donald Curtis started the project.

For the historians in the audience, Donald has previously worked with me on Ghemical-gms. Ghemical is an open-source model builder based on Babel written by Tommi Hassinen, which Donald made installable on Windows and Mac's to which he added an interface to GAMESS. During the process Donald and Geoff hooked up and the rest, as they say, is history.

Sunday, October 18, 2009

Proteopedia is just what the name suggests: a wikipedia for proteins. What makes it specific to proteins is the integration of Jmol for animated and interactive figures. It has a long list of pre-made entries ready for use in a lecture and, like a real wiki, you can add content as well.

A lot of the features are explained using a set of very nice screencast. The screencast above is simply a recording of the first 2 minutes of the first video (protopedia does not offer a way to embed the videos) narrated by Eran Hodis.

Wednesday, October 14, 2009

For a while I thought I would be only science blog that would make it through October without mentioning the Nobel Prize. But via Chemical Blogspace I came across this video interview of the 2009 Nobel Prize co-winner in Chemistry Venki Ramakrishnan posted on The Sceptical Chymist (note: spell-check when you name your blog!):

During the interview they show some stunning molecular animations of the ribosome in action, which are credited to Said Sannuga and Ramakrishnan. I wonder ... Sure enough there are several very interesting animations on Ramakrishnan's web site. The screencast above shows a few seconds of one of them. For the rest, go to the site and have a look.

Sunday, October 4, 2009

In a previous post I showed how to use PubChem Editor to draw a complicated molecule and import it into Avogadro using the insert SMILES option. Unfortunately, Avogadro does not interpret the chirality information correctly, and chiral structures often have to be fixed.

In the screencast above I show another program, Marvin Sketch, that allows you to draw a structure in 2D and convert it (correctly) to 3D. The coordinates can then be transferred to Avogadro by saving a .mol file.

I also use Marvin to draw 2D structures for the blog. For example

However, Marvin Sketch is much more than a drawing program. The tools menu contains many powerful cheminformatics tools such as pKa prediction, isomer search, etc.