In this tutorial we will be optimizing the reactants, products, and the transition state of a reaction. We will learn how to prepare Gaussian 16 input files with the objective to get its potential energy surface.

DFT Calculation Workflow for Optimization

Before tackling the transition structure (TS), the first step is to optimize the geometries of reactants and products. If the reactants or products are not stable minima, then it is impossible to locate a TS that connects them.

1. Build the molecules using Gaussview (or use XYZ coordinates from another source) to create an input file for Gaussian (.com file). Below is an example output file for a water molecule

%mem=8GB
%nprocshared=8
%nosave
#P METHOD/BASIS-SET opt freq=noraman

opt freq=noraman water

0  1
O          0.24600        0.20600        0.00000
H          1.20600        0.20600        0.00000
H         -0.07500        1.11100        0.00000

General considerations:

A Gaussian input file has a specific format which must be followed. Blank lines are not optional and must be there at specific positions. Actual example input files are provided below. Some options can be included to specify how the calculation will be run on the computer/cluster:

%mem=XX tells Gaussian how much memory to use
%nprocshared=XX tells Gaussian how many shared processors to employ
%nosave tells Gaussian not to save the intermediate files (checkpoint .chk, integrals .int, read-write file .rwf, etc.) if the calculation is successful. For some calculations it might be required to save some of those files and to change this option.

#P starts the route line, where you indicate the job you want to run. Aspects needed:
METHOD: The computational method (for us usually density functional) used
BASIS-SET: The basis set used
EXAMPLE: # B3LYP/6-31+G(d,p)

JOB COMMANDS: Keywords to tell Gaussian what we want to do.
For ground-state optimizations (like reagents, products, intermediates): opt freq=noraman 
For TS optimization: opt=(ts,calcfc,noeigentest) freq=noraman 
–opt requests a structure optimization, freq=noraman requests a frequency calculation on the optimized structure. While you can split these two calculations in two input files, we find it easier to have them together in one output for anything but the largest calculations.–

After the route line, you need a blank line, then a title line (it can be anything), then another blank line.

Then, you place the information about your structure. First, you need to specify the charge of the system (0 for neutral, -1 of monoanionic, +1 for monocationic, etc.) and the multiplicity (typically 1, for singlet). Finally, you insert the XYZ coordinates of the structure, followed by a final blank line.

2. Transfer the file on your local cluster with a file transfer software (we use FileZilla).
3. Submit the input file to a Gaussian calculation (how to do this varies for each cluster)
4. Once the job is finished, download the output file (.log or .out) and visualize the results using Gaussview or another visualization software.
5. You need to verify that the structure has no imaginary frequencies, so it can be considered a true minimum. Then, you can extract the energies, enthalpies and free energies from these output files.

Workflow for Transition State Localization

A transition state is the maximum of energy along the reaction coordinates (seen in top graph), and a minimum in all other directions (1st order saddle point) (seen in graph below). Locating such structures is both a science and an art, and various approaches are possible. Two main approaches are 1) to run a scan along the reaction coordinate in order to find a structure close to a maximum in energy or 2) to run a constrained optimization on a guess structure. Both approaches require running a constrained optimization in Gaussian.

Constrained optimizations to find transition structures:

• In order to locate the TS using a scan, you want to progressively modify the structure of a minimum connecting the TS (a reagents complex or a products complex) along the reaction coordinates and calculate the energy of each structures to find a maximum.
• If you cannot run a scan and do not have access to a good guess structure from another source, you will have to optimize a guess geometry before doing the full TS optimization.

Gaussian can optimize a structure while keeping a geometry parameter (bond length, angle, dihedral) frozen. Bond lengths/angles that are critical for the reaction coordinate (e.g. a forming or breaking bond) should be frozen so that Gaussian does not try to optimize them (i.e. make them at the length of a formed/broken bond). Alternatively, you can request Gaussian to scan along that geometry parameter that is critical to the reaction. Below is an example of a scan input file.

%mem=2GB
%nprocshared=8
%nosave
#p opt=modredundant b3lyp/6-31+g(d,p)

scan SN2 F-

-1 1
 C                  0.00000000    0.29534700    0.00000000
 Cl                 1.84578500    0.28507800    0.00000000
 H                 -0.30677300   -0.27161000    0.87777800
 H                 -0.30677300   -0.27161000   -0.87777800
 F                 -2.85029100   -2.51104600    0.00000000
 C                 -0.51958200    1.71858400    0.00000000
 H                 -0.18959400    2.26610200   -0.88878200
 H                 -1.61550300    1.68053000    0.00000000
 H                 -0.18959400    2.26610200    0.88878200

B 1 5 S 35 -0.1
A 5 1 6 F

The format of the input file is similar. You request a constrained optimization by requesting the job type opt=modredundant. After the blank line below the XYZ coordinates, you have to specify your constraints, i.e. what parameters will be fixed or scanned.

To constrain:
B 1 2 F: freezes the bond length between atoms 1 and 2 (“B” for bond)
A 1 2 3 F: freezes the angle between atoms 1, 2 and 3 (“A” for angle)
D 1 2 3 4 F: freezes the dihedral angle between atoms 1, 2, 3 and 4 (Z) (“D” for dihedral)
N.B.  Gaussview displays the atom numbers when clicking on them 

To scan, write your constrain as: B x y S N s
“B”: bond (it could also be an angle or dihedral, you simply need to specify more atoms as needed)
“x y”: between atoms #x and #y 
“S”: scan 
“N”: number of steps for the scan
“s”: step size (in Å for bonds or degrees for angles) 

In the above example, you can notice we scan a bond between atoms 1 and 5 for 35 steps of -0.1 Å. We also keep the angle between atoms 5, 1 and 6 fixed.

Once the scan is complete, the structure corresponding to a maximum of energy (and low gradient) should be very close from the actual TS for this system. This structure can then be submitted to a full TS optimization (as explained above).

For transition structure, then should have exactly one imaginary frequency to be considered true first-order saddle points.

Examples

As an example, let’s compare the S_N2 and E2 reactions of chloroethane with fluoride.

Optimizing the reactants and products:

Reactants:
Chloroethane: (.com input file)
Fluoride: (.com input file)

S_N2 Products:
Fluoroethane: (.com input file)
Chloride: (.com input file)

E2 Products:
Ethylene: (.com input file)
HF: (.com input file)
Chloride: (.com input file)

Optimizing transition structures:

TS for S_N2: (.com input file)
TS for E2: (.com input file)

For each structure, the input file is provided, so you get to see examples of each type. In case you did not have good guesses for the transition states, you would want to locate them yourself.

Creating a guess structure:

• Prepare a guess structure where both the breaking C-Cl bond and the forming F-C bond are 2.2 Å. Then, ask Gaussian to optimize by keeping those two bonds at this distance. This choice is arbitrary but most be informed by chemical intuition or similar structures.
• In this specific case, atom 1 is the carbon atom linked to the chloride in chloroethane, 2 is chlorine, 5 is fluorine.
• Two constraint lines at the end of the file are to freeze the C-Cl bond (B 1 2 F) and the C-F bond (B 1 5 F) at the value they have in the input file. 
• After the optimization is finished, the final structure can be submitted to a full TS optimization.

Constrained optimization of an S_N2 guess structure: (.com input file)

Generating a scan:

• The reaction coordinate of an S_N2 reaction involves simultaneous bond breaking and bond forming. You have to choose one bond to scan, so let’s choose the F…C distance.
• Starting from the reactant complex, we will scan the bond between carbon 1 and fluorine 5 by increments of 0.1 Å.
• The input file also freezes the angle between fluorine, carbon and chlorine.

Scan along the S_N2 reaction coordinate: (.com input file)

Results from the scan:

• The scan starts with the fluorine too far from the carbon, so the first points decrease the energy. Around point #12 is the reactant complex, and from there the energy rises to a TS at around point #19. The product complex is around point #27.
• Structure #19 is the guess structure for the TS, which will be submitted to a full optimization.

Final results:

Some Issues and Their Solutions 

For TS optimizations:

If during a TS search the structure optimizes toward a local minimum, with an important decrease of the energy on the optimization graph, that means the starting structure for the TS search was too far from the actual TS geometry (and consequently too close from a minimum). The solution is to start from the submitted structure (first point) and modify the structure to make it further from the minimum. For example, if you optimize a TS for a bond formation, and the TS search breaks into a structure where the bond is formed, restart the optimization from the beginning with a longer distance between the two bond-forming atoms.

For any kind of optimization:

Often you get a Gaussian error called “convergence failure”. To navigate around those, two main solutions have worked for us:

• solution #1: from the lowest energy point, slightly modify the structure (change a angle or a bond length) and reoptimize.
• solution #2: from the lowest energy point, reoptimize the structure with a different SCF (self consistent field) method (keyword: SCF=option). Usually, SCF=QC is useful for difficult convergence cases.