1. No category

Transition State Searching in Gaussian/GaussView – an SN2 example
This tutorial is an example of transition state searching. Gaussian does the actual calculations while
GaussView is the build/edit/analysis tool. Consult our "Introduction to Quantum Mechanics – Basic
Calculations using Gaussian" tutorial for basic molecule building and calculation setup instructions.
Build reactant and product
Start GaussView and click on the 6C icon. Select Si in its trigonal bipyramidal geometry (the last
geometry on the right), and click in the black window (workspace). Now select C in its "C Atom"
geometry and click on the Si atom in the workspace. Select Cl and click on one axial hydrogen, then do
the same using Br for the other axial hydrogen. Select the Modify Bond icon, then in the workspace
click on C then Cl. In the dialog that comes up, set the Atom 1 pulldown to Fixed,
and set the bond length to 1.8 (see Figure 1). Do the same to the C-Br bond but set
it at 2.4. Select the Modify Angle icon, and click on Br then C then an H atom.
Again, set Atom 1 to Fixed (Atom 2 is also Fixed, Atom 3 is Rotate Group) and
set the angle to 70.0. Repeat this step to set all three Br-C-H angles to 70.0. You
must click OK each time to make the changes take effect. Otherwise they will reset
when you close the Modify Angle window. You have just built an approximate geometry of the reactant
but you need to optimize it fully in Gaussian to obtain reliable energy, MO, and related data.
Select Calculate - Gaussian Calculation Setup and select Opt+Freq from the first pulldown in the Job
Type tab. In the Method tab, make sure the Method is Ground State - DFT - Default Spin - B3LYP.
Set Basis Set: to LanL2DZ. Set the charge to -1 and press the Enter key. The Spin should
automatically change to Singlet. The charge and multiplicty must be correct or the calculation is
meaningless! Click Submit and save the file as greact.com when prompted. Use the qstat command to
check on job status. If your job is not listed, it has either died or finished successfully.
While that greact calculation is running you can set up the product structure. Select File - New - Add to
Molecule Group. Your workspace now has two views (or entries), the first with your reactant and the
second which is empty (toggle between views using the arrow buttons at the top of the black window).
With your reactant displayed in the first view, select Edit - Copy, which places
the structure in the Builder Fragment window. Toggle to the empty second view
and click in the black space to place a copy of the reactant there. Adjust the bonds
and angles in the second view to invert to reactant into the product. Do not just
change the Cl into Br and vice versa – the atom ordering is important for vector
following in transition state calculations so you need to modify reactant into
product by adjusting existing bond lengths and angles. The C-Cl and C-Br distances should be reversed
(see Figure 2) and the Cl-C-H angles are 70.0 degrees. With this new structure (product) displayed, set
up and submit a calculation the same way as before, saving the file as gprod.com.
Viewing optimization results
Once qstat shows that your greact job is no longer running, open greact.out using a text editor like gedit.
At the end of the file you should see "Normal termination of Gaussian 09 at ..." (if the job finished
successfully.) Search in the file for the phrase: Optimization completed. Check the 4 convergence values
and thresholds in the lines just above. Next, search for the phrase: Harmonic frequencies. This section has
the vibrational modes. Note that the first few Frequencies are all positive indicating the structure is a
minimum on the potential energy surface. Scroll down to the thermochemistry section and find the line
with the phrase: Sum of electronic and thermal Free Energies. Write down the corresponding number
including all of the digits to the right of the decimal and also note the units for that number.
In GaussView, close all existing displayed molecule windows and then go to the File - Open dialog.
Check the box: Read Intermediate Geometries at the bottom of the dialog, then open greact.out. This
loads in the initial geometry and all of the partially optimized geometries found as part of the
optimization procedure. The default displayed structure is the initial input! Use the arrow buttons at the
top of the workspace to click through the intermediate steps till you get to the final optimized structure.
Does this look like a reasonable geometry for the reactant? Then, in the Results menu, select
Optimization to plot a graph of electronic energy versus optimization step number. You can also navigate
through the geometries by clicking on graph points. To zoom in on an area of the graph, drag out a
rectangle with the mouse. How does Gaussian know when an optimization is complete (converged)?
[Geometry convergence failures are common in QM in that the geometry may be fine but fails the built in
convergence test. This requires manual intervention. You would evaluate the final geometry to decide if it
is reasonable and also visually inspect the optimization graph to see if the energy is mostly converged.
Then you would resubmit your Opt+Freq as a Frequency only, starting from that geometry.]
With the last structure displayed in the black window, select Results - Vibrations. This displays the
vibrational frequencies and their IR intensities. Visualize some vibrations by selecting a vibration and
clicking the Start Animation button. [Click Stop Animation before closing the Vibrations dialog or it
keeps going.] Next, do the same steps on gprod.out in GaussView and gedit and analyze those results.
What is the G of this reaction? Note: for meaningful energy comparisons, structures must be stationary
points: minima with real (positive) vibrational frequencies, or TS with one imaginary (negative)
frequency. [If you don't have frequency data (necessary for H and G), you can compare the electronic
energies of structures (listed as E(RB3LYP) in Results - Summary. This number can also be found in
greact.out. Look for the last instance of a line labelled: SCF Done.]
Setting up and running the transition state search
In GaussView, close all displayed molecule windows and then go to the File - Open dialog. Un-check
Read Intermediate Geometries and open greact.out. Go back to that dialog and change the Target:
pulldown to Add all files to active molecule group and then open gprod.out. The two structures must
be in one molecule group (two views, not two separate windows). Toggle between the views to verify that
one is reactant and one is product. Select View - Labels and make sure the atom numbering is identical
in reactant and product. Click the Modify Bond icon in the GaussView mainframe, then click the C and
Br atoms and select the bond order between C and Br as half-bond (dotted line) in the pop-up dialog box.
Repeat this and set all C-Cl and C-Br bond orders as half-bond in both the reactant and product. Open the
Gaussian Calculation Setup dialog and change the Optimize to a: pulldown to TS (QST2) (Job Type
should still be Opt+Freq). The Keywords area (top of the dialog, also called the Route section) may list
a bunch of options we don't want, this is GaussView trying unsuccessfully to help out. Delete anything in
the Additional Keywords: field on the Job Type tab, set Guess Method: to Default on the Guess tab,
and finally set Model: to None in the Solvation tab. The Keywords should ONLY contain "# opt=qst2
freq rb3lyp/lanl2dz geom=connectivity". Submit the calculation as gts.com.1
When the job finishes (a few minutes), close all displayed molecule windows and open gts.out (with
"Read Intermediate Geometries" set). Inspect the vibrational frequencies to verify you found the correct
transition state: there should be one imaginary ("negative") frequency which corresponds to the reaction
coordinate (what is the meaning of a result with more than one imaginary frequency?). Compare the free
energy of the transition state with the reactant and product. Is the TS reactant-like or product-like?
The bond types shown in GaussView do not affect the Gaussian calculation. However, for some structures,
GaussView cannot read a checkpoint file (.chk) without a “geom=connectivity” keyword and bonds between
atoms with long distances and weak interactions. This is a bug in the current version of GaussView.
1
It is often clear from inspecting the transition state's vibrational frequencies whether it corresponds to the
correct reaction coordinate, but for publication quality calculations it is necessary to prove this by
performing an Intrinsic Reaction Coordinate (IRC) scan (see below), which links the transition state
structure to those of reactant and product.
Relaxed coordinate scan
While the above calculations completed properly, sometimes finding an initial guess for a transition state
is more difficult. If you think the reaction coordinate is mainly composed of a bond stretch or angle bend,
a relaxed coordinate scan is a good way to approximate the potential energy profile by generating a series
of structures along that path. This method performs a geometry optimization of all coordinates except one
that is specified, then changes the specified coordinate and optimizes again, etc. as many times as you
specify by giving the number of steps and an increment size.
*If you are doing this as a class exercise, you can skip running the coordinate scan and IRC. The
following instructions are mostly here for when you need them for your research. For now, go to
the last page and do the section on displaying electrostatic potential. Try also to export/import
structures from other programs as described on that page, and also look at the web links.*
Open gprod.out in GaussView and make sure you have the final structure displayed. Select the Modify
Bond icon, fix the C atom and modify the Br-C distance as 2.1 Å. Right click in the (black) workspace
and select View - Labels. Select Redundant Coordinates from the Edit menu. In the Redundant
Coordinate Editor, click the Add button (near the top), and change the pulldown listing "unidentified" to
Bond. In the two Coordinate: fields type in the atom numbers for the Br atom and the C atom. In the
pulldown listing "Add," select Scan Coordinate. Now, fill out the fields to the right so that the phrase
reads "Take 24 Step(s) of size 0.05 angstroms." This will result in a C-Br distance scan from 2.1 Å to 3.3
Å in increments of 0.05 Å. Keep the field listing "Don't Set" as it is (DO NOT reset the starting value of
the scanned coordinate. If you need to start the scan with a different value, modify the structure.). Click
Ok. In the Gaussian Calculation Setup dialog, verify that the pulldowns in Job Type are set as Scan and
Relaxed (Redundant Coord). Take out any inappropriate keywords so that your Route section looks like
this: “# opt=modredundant rb3lyp/lanl2dz geom=connectivity”. Now submit the job (either from
GaussView or from the command line), calling it gscan.com.
When the job finishes, open gscan.out in GaussView with Read Intermediate Geometries selected.
There should be 25 structures total. Select Results - Scan and inspect the energy profile graph.
Intrinsic reaction coordinate (IRC) scan
The IRC scan needs to know the force constants at the start of the job. Although you can read the force
constants from a frequency calculation in your IRC calculation, it is more reliable to calculate the force
constants at the beginning of an IRC calculation.
In GaussView, open gts.out and make sure you have the final structure displayed. Select Calculate Gaussian Calculation Setup. In the Job Type tab, select IRC from the first pulldown, and then set the
Force Constants pulldown to Calculate Once. Notice that Follow IRC is set to Both directions, which
means it will go both ways from the TS along the reaction coordinate. Select Compute more points, and
change the N= field to 30. This is the number of points along the reaction path to examine (in each
direction if both are being considered). Thus the calculation in each direction will stop after finding 30
points if it hasn't already found a minimum. Delete inappropriate keywords so that your keywords are
only: “# irc=(maxpoints=30, calcfc) rb3lyp/lanl2dz geom=connectivity”. Submit this as girc.com.
Monitor the job using qstat as usual. In GaussView, you can view the structures in the reaction pathway
by opening girc.out (with the Read Intermediate Geometries option) and selecting Results - IRC/Path
and inspect the energy profile graph that is displayed. At the end of output file, you will see “This type of
calculation cannot be archived”. This is common for some normal terminated IRC calculation.
------------------------------------------------------------------------------------------------------------------------------Other Common Gaussian Tasks:
Instead of submitting calculations directly from GaussView, you can click Retain at the bottom of the
setup dialog, then go to File - Save in the main GaussView window and save a file as filename.com. You
can open filename.com with gedit to inspect keywords or make other adjustment. You can submit the job
by typing on the terminal command line (make sure you cd to the directory containing filename.com):
run_gaussian filename.com
Displaying electrostatic potential surfaces
It can be useful to display the charge distribution projected on a molecular surface. To do this, open
filename.chk.and make sure you have the final structure displayed (if multiple structures were read in).
Select Results - Surfaces/Contours, and click on Cube Actions - New Cube. Change the Type:
pulldown to Total Density and click Ok. In a couple of seconds (longer for a larger molecule) an entry
will appear in the Cubes Available section of the Surfaces and Contours window. Repeat this step but set
Type: to ESP. Now, select the Electron Density from Total SCF Density entry in the Cubes Available
section of the Surfaces and Contours window, change the Density= field to 0.001, and then click on the
Surface Actions button and select New Mapped Surface. Click Ok in new dialog. You should now see a
surface obscuring your molecule in shades of orange. On the top left field above the workspace window,
change the value to -0.2. Now, right click in the black space and select View - Display Format Surface and change Format: to Transparent. Adjust the transparency level using the slider.
Exporting structures from Maestro to GaussView (and back again)
Maestro to Gaussview: In your Maestro Project Table, make sure the structure you want to export is
highlighted in yellow, because by default the export tool chooses selected entries, not what is in the
workspace. Select Project - Export Structures. Set the Files of type: pulldown to MOL2 (this format
usually works but PDB is a good alternative choice if the mol2 doesn't work well). Name the file in the
File name: field and click Save. Open GaussView and select File - Open. Change File type to MOL2
Files, select the file you want and click Open, and then click Ok in the new dialog.
Gaussview to Maestro: In GaussView, select File - Save and change File type: to MOL2 (.mol2) or PDB
(.pdb). Be sure to type the file suffix along with the filename in the File name: field. Click Save. In
Maestro, use Import Structures to read in the structure. If the mol2 or pdb file you created in
GaussView cannot be imported, try opening it in Mercury (type mercury on the command line) and recpnvert the file there. Mercury is also the best choice to display X-ray crystal structure files (.cif).
Links
Please see our webpage and the documentation for more advanced Gaussian09 tasks:
http://glab.cchem.berkeley.edu/glab/chemmisc.html#gaussian
Advanced input files, restarting, etc:
http://glab.cchem.berkeley.edu/glab/faqs/gaussian_part2.html
How to understand thermochemical data in Gaussian:
http://www.gaussian.com/g_whitepap/thermo.htm