How to simulate chemical reactions in GAMESS?

I am rather a newb in the GAMESS field but finally, learned to make input files for most molecular simulations and could even use TD-DFT for excited states (special thanks to Geoff Hutchison for his help).

Now I want to start simulating simple chemical reactions. I looked at the GAMESS manual and think I need to make transition states and find the least energy pathways but have no idea how to make an input file for that.

Can anyone give me a simple input file example for a simple reaction?

As in geometry optimization, you are searching for a stationary point on the potential energy surface (PES). Not for local minimum, but for saddle point, therefore in GAMESS, you specify RUNTYP=SADPOINT. You would also need the correct (non-guess) Hessian matrix, which you can calculate separately with RUNTYP=HESSIAN. In the transition state (TS), you should have zero gradients (stationary point) and one imaginary eigenvalue of Hessian (the reaction coordinate in TS).

For a working example, see e.g. SN2 transition state. In principle, you need a very good guess at TS geometry, which you obtain by scanning along some chosen coordinate. The point highest in energy should be reasonably similar to the TS, and you optimize to a saddle point starting from there. Already at the guess geometry, you must have one imaginary Hessian eigenvalue (resembling the reaction you are looking for, not methyl rotation somewhere else).

But beware, finding the transition states is one of the most difficult topics in computational chemistry and you are not guaranteed to find one, no matter how hard you try. One of the biggest problems is solvation, as you often wish to find a reaction barrier in solvent (where the majority of the chemistry happens), but you can calculate only in a vacuum. Just think of how carefully the workbench chemist finds the correct solvent for a given reaction, on the computer you have only very crude approximations to it (and most often just vacuum).

That said, it is great fun and rewarding to find the transition states. For best results, use the simplest theory available, as it is better to have smooth and well-behaved PES, even though it is slightly wrong, so it is no harm to start with HF or BLYP (or B3LYP). To get some ideas, start with simple reactions, i.e. where no charged species are involved, so electrocyclic reactions are probably the best choice, e.g. Diels–Alder.

• (+1) I am not an expert in GAMESS, but what you say makes sense in general. As for the procedure, I agree with what you say. I would recommend to start with BP86, because it is a pure functional, quite robust and pretty fast. Best also use a basis set of double zeta quality, like cc-pVDZ or def2-SVP. Better avoid the quite old Pople basis sets, like 6-31G*. Also use a fine integration grid and tight convergence criteria. It takes longer, but usually the optimizations don't run in a undesired direction. A Diels-Alder reaction is certainly a good training exercise. Commented Nov 28, 2014 at 10:36
• Thank you ! This was a great answer and I learned a lot with that. As I realized from that SN2 Transition State example, we should already know ( or at least guess ) the result and then make our energy path on that. Is there any way to do it in reverse? I mean doing the simulation to make some crude guess about what "may" happen in a new reaction ?
– Aug
Commented Dec 4, 2014 at 0:05
• @Aug: in general, you cannot ask computer: "I have those two molecules, how will they react?" If you assume that new bond will be formed between them, than the computer should try to propose MxN reaction, where M, N is number of atoms in each molecule. But what if it is electrocyclic reaction? Or elimination? So it is up to you, chemist with brain, to propose something and ask computer: "Does this reaction look reasonable?" and get answer "Yes/No/Depends" Commented Dec 4, 2014 at 6:44
• People are starting to put the results of their calculations on figshare.com, eg chemapps.stolaf.edu/jmol/jmol.php?source=http://… To see the reaction, (i) right click somewhere next to the molecule and a menu will pop-up; (ii0 go to animation -> animation mode and press loop; (iii) right click somewhere next to the molecule to get the pop-up menu and press play. There are many other examples. Commented Oct 12, 2015 at 12:06