Gaussian: Common Topics and Useful Procedures

This file picks up after our Gaussian Intro tutorials. We have some other Gaussian related help files so if you don't see your answer below, please also check the faqs page, our notes on custom basis sets, our notes on convergence problems, and our main software page. See the Gaussian manual for more information.

This provides information about commonly-used operations when performing quantum chemical calculations in Gaussian. The information in this how-to, together with that in the basic tutorial, should provide you with the BASIC knowledge necessary to answer normally-occurring chemical questions. Some information regarding theory will be applicable to other QM programs. The later sections are explanations of operations not fully addressed in the program documentation. Feel free to make suggestions for additions and/or corrections to this file.

Contents:

Creating and Submitting Input Files, Advanced Version

Starting Jobs from Checkpoint Files

Vibrational Frequencies

Transition State Searches

ONIOM Calculations

G3 Energy Calculations Using Previously-Calculation Structures

NBO (Natural Bond Orbital) Analysis

Obtaining Graphics Files

High-Quality MO (or NBO) Pictures


Creating and Submitting Input Files, Advanced Version:

While many of the commonly-used options are available from Gaussview menus, you may find that you need to manually add something to (or delete something from) the input file. You can also construct new input files without ever entering Gaussview.


Saving files:

To save an input file in Gaussview without submitting the job, click on Retain in the Gaussian Calculations Setup menu instead of Submit. You can then save the file through File - Save.


Submitting from the command line:

From the directory in which your input file is located, submit the job by typing run_gaussian filename (with or without the .com extension) in the command line. The job will now be running in the same fashion as if submitted via Gaussview.


Input file format:

The Gaussian input file is text-based and easy to edit manually. A basic file is composed of four parts, the Link0 section (lines starting with %), the Route section (line starting with #), the job title, and the Molecule section. The placement of blank lines is important, so follow the example below if you construct an input file from scratch. Following is a sample input file for a DFT (B3LYP) geometry optimization on water, also requesting a vibrational frequency calculation, using the 6-31G(d) basis set:

%chk=/home/username/directory/water1.chk
%mem=1GB
%nproc=1
#T opt freq rb3lyp/6-31g(d)
{blank line}
water optimization
{blank line}
0 1
O
H                  1              B1
H                  1              B2  2              A1

B1             1.37000000
B2             0.96000000
A1           109.50000006
{blank line}

The Link0 section specifies options related to computer usage and file locations. The most important line here is the %chk= line which specifies the name and location of the checkpoint file that the calculation will create. If only the filename is given (e.g. %chk=water1.chk), it is assumed it is located in the directory from which the job was submitted. The other two lines relate to processor and memory usage, and should not normally be altered unless you already know the details of how the queuing system works (again, refer to the Faqs for using more than one processor).

The Route section specifies non-default computational details, which can be written in any order. For instance, if we put freq before opt in the above file, the optimization would still proceed before the frequency calculation. More than one line can be used, but the section must always terminate with a blank line. The T after the # character is optional and specifies terse output (or P for verbose output). This does not affect what is calculated and generally aides in locating important information in filename.out. Many keywords have arguments associated with them; for example, if we wanted to make the convergence criteria for the optimization more stringent we could enter opt=tight or opt(tight). For multiple arguments, use the syntax opt(tight,maxcycle=200) which would also tell the program to run a maximum of 200 optimization iterations. The list of possible keywords and arguments is large, so refer to the Gaussian manual for more information.

The first line of the Molecule section specifies only the charge and multiplicity (2S+1) of the molecule. Here we have 0 1 for neutral singlet water. In this example the coordinates of the molecule are listed as a Z-matrix, which is very useful if you want to constrain specific coordinates (follow the coordinate with a space and then the F character). You could also use Cartesian (x y z) coordinates.


Starting Jobs from Checkpoint Files:

It is often desirable to restart a failed calculation (see separate help file) or start a new calculation using the structure, wavefunction, and/or the force constant matrix from a previous calculation. For example, if SCF convergence was a problem you might run a single point calculation using keywords, meant to help SCF convergence, that would hurt efficiency in a geometry optimization. You would then start a geometry optimization specifying that the previously-determined wavefunction should be read from the checkpoint file as an initial guess.
Note: if you change the geometry of the structure or add/delete any atoms, you cannot use previously-calculated wavefunctions or force constant matrices from a checkpoint file! However, you can add/subtract electrons, just remember to correspondingly change the multiplicity.

Start by making a copy of the checkpoint file (cp original.chk restart.chk) of the failed/failing/previous calculation. The copied checkpoint filename should reflect the name of the new calculation you are about to start from the failed/failing/previous calculation state. This prevents overwriting the old checkpoint file (which you may need later on) while allowing the restart to get needed information from the previous calculation. Read the copied checkpoint file into Gaussview and start where you left off by saving a new input file (restart.com) matching the new checkpoint file name.

Some options for each keyword that specify reading from the checkpoint file are: Opt=ReadFC (reads force constant matrix if it existed in the checkpoint file, also works with Freq, see below); Guess=Read (reads the wavefunction from the checkpoint file, uses it as an initial guess for SCF); Geom=Checkpoint (reads the structure from the checkpoint file, you still need to specify charge and multiplicity); Geom=AllCheck (reads the structure, charge and multiplicity from the checkpoint file); Geom(Step=N) (reads the Nth structure from a previous geometry optimization checkpoint file). The last three arguments are useful if you are constructing the new input file in a text editor, as it eliminates the need to copy/paste atomic coordinates.


Vibrational Frequencies:

If you have run an optimization (minimization or transition state search) you can request a vibrational frequencies calculation (keyword Freq). Even if you do not need the values of the vibrational frequencies or thermodynamic data, it is still necessary to verify that the structure is a stationary point (dE/dr = 0 where r is the collection of all coordinates) and is of the correct type. This can be determined by the following: 1) a minimum is characterized by no imaginary frequencies (listed as negative frequencies in the output); 2) a transition state is characterized by EXACTLY one imaginary frequency, which should correspond to the reaction coordinate. More than one imaginary frequency indicates a higher-order saddle point and normally has no physical meaning.
Note: imaginary frequencies of small magnitude (less than 10 wavenumbers, or larger according to some) can usually be ignored.

You can specify a frequency calculation along with an optimization, as in the example above. The frequency calculation will automatically proceed when (or if) the optimization completes. Alternatively, you could start a frequency calculation using the checkpoint file from a previous optimization. In this case it is a good idea to use the keyword Guess=Read, so the wavefunction does not need to be recomputed. You can specify Freq=ReadFC so the calculation uses the approximate force constant matrix constructed during the optimization, however the analysis (especially the zero-point energy) will be of poor quality. It can be used as a first-glance check because it is a very fast calculation, but you should go back and recalculate the frequencies analytically (the default for HF and DFT methods) before placing any confidence in the results. You can view and animate the vibrational frequencies in Gaussview by reading in the output (.out) file from the frequency calculation and selecting Vibrations from the Results menu.

In addition to the vibrational frequencies, this calculation will also calculate IR intensities, and molecular thermodynamic properties (at whatever temperature/pressure you specify, default is STP). The output lists corrections to the electronic energy including zero-point energy, thermal energy, enthalpy, entropy, and Gibbs free energy, all in hartree (1 hartree = 627.5095 kcal/mol).

The choice of which of these corrections to include is currently still a matter of debate. Obviously the Gibbs free energy should most closely reflect experimental findings, yet errors associated with calculating the vibrational frequencies can make these figures somewhat unreliable. A common compromise is to use electronic energies corrected only for zero-point energy, which can be viewed as 0 K enthalpies.


Transition State Searches:

In a broad sense, a transition state search is the same as a geometry optimization except that one coordinate is maximized instead of minimized. Here the term coordinate actually refers to an eigenvector of the Hessian (force constant matrix), which can be a combination of bond lengths, angles and torsions. This eigenvector corresponds to the reaction coordinate and has a negative eigenvalue in the Hessian. It is this negative eigenvalue that gives rise to the one imaginary frequency obtained in a frequency calculation on a transition state (you should now have a feel for the meaning of a frequency of imaginary amplitude).

There are three methods of setting up a transition state search: TS(Berny) (keyword Opt=TS), TS(QST2) (keyword Opt=QST2), and TS(QST3) (keyword Opt=QST3). The TS(Berny), or standard, search only requires that you input a guess of the TS structure. You must have all or part of the Hessian calculated and supply this as an initial guess, as a standard initial-guess Hessian is insufficient to identify the reaction coordinate. Obviously, your initial guess structure must be a very good guess and so this method is generally a bad choice when starting a calculation from scratch. It is useful, however, for restarting a calculation using a structure that appeared "on the way" to optimizing correctly but then took a turn for the worse.

Setting up the standard TS search in Gaussview is nearly the same as for normal geometry optimization, selecting Job Type -- Optimization -- Optimize to a: TS(Berny). You should also set Calculate Force Constants: Once (Opt=CalcFC). This calculates a good Hessian in the beginning so the reaction coordinate can be located.

The QST2 method (similar to LST in other programs and Gaussian98) is generally the most useful of the three ways to start a TS search. Here you input two structures (see below), one of the reactant and one of the product. The program then interpolates between this structures to get an initial TS guess, and information about the interpolation is used to help identify the correct Hessian eigenvector as the reaction coordinate. In practice, it is usually best to help the program by inputting reactant and product structures intermediate between the minima and the expected transition state. This tends to isolate the major components of the reaction coordinate near the transition state and improve the chances of correctly identifying it.

The QST2 method setup is more complicated because you need to enter two structures, which must have identical atom numbering schemes. You do this by using the Molecule Group (Molgroup) function. The order in which you specify reactant and product is irrelevant due to the principle of microscopic reversibility. The easiest way of obtaining both structures with identical atom numbering schemes is the following: 1) draw or import the reactant (or product); 2) select Copy from the Edit menu; 3) select Paste -- Add to Molecule Group (you now have two identical structures which you can toggle between using the top left of the molecule view window); 4) edit the second structure using bond length and angle adjustments (you may need to make and/or break some bonds) so it corresponds with the product (or reactant) structure; 5) select Job Type -- Optimization -- Optimize to a: TS(QST2). If the TS search fails by falling to one of the minima, say the product, you need to go back and edit the product structure so that it is closer to the expected TS structure. This effectively shifts the interpolated TS guess closer to the reactant.

The third method, QST3, is a combination of the previous two, wherein you specify reactant, product, and a TS guess. This can be helpful in guiding the reactant-product interpolation in cases where the QST2 method fails to give an adequate guess. The QST3 method is set up the same as the QST2 except that you specify a third molecule, which is the TS guess and must be molecule #3 in the Molecule Group. Most of the time this TS guess will be a previous attempt that is "almost there." It can be added to the molecule group by doing the following: 1) read in the .com file for the QST2 job used to get the TS guess; 2) make sure structure 2 is the structure selected and displayed on-screen; 3) import the TS guess-containing file, after setting the Target field in the Open Files dialog to "Add all files to active molecule group."


ONIOM Calculations:

The ONIOM method allows the use of more than one model chemistry within a molecule or system of molecules. It is useful in large systems where a high level of theory is needed for one or more regions and the remaining portion can be treated at a lower level. This is often a better alternative than using a simplified system, as the latter often involves sacrificing important steric features. The levels of theory can be any mixture of QM and/or MM models. See Dapprich, J. Mol. Struct. (Theochem) 1999, 461, 1 for an introduction to ONIOM and some important considerations. Read that first, as the following discussion uses terms defined therein. Gaussian allows for the use of two or three layers of theory, but this how-to only covers the two-level treatment. ONIOM calculation are best set up in Gaussview unless you really know what you are doing.

Setting up the layers:

Build or import your molecule in Gaussview. Select Select Layer from the Edit menu, which brings up the Layer Selection Tool. With Low in the Set Layer field, click on the Select All button and then Apply. The molecular representation will switch to wireframe. Select Set Layer -- High. While holding the Shift key, click on each atom you want included in the High layer (higher level of theory), then Apply. The High level atoms should switch back to ball and stick.

Setting up the calculation:

This is done the same way as normal except that you must check the box next to Multilayer ONIOM Model in the Method tab of the Gaussian Calculation Setup dialog. There should now be three sub-tabs, one for each layer (layers that were not specified are greyed out). Specify the method details for each layer.
Note: for a two-layer job, there will be three entries for charge/multiplicity. The first set is for the complete system (Reallow), the second is for the model system at the high level of theory (Modelhigh) and the third is for the model system at the low level of theory (Modellow). Modelhigh and Modellow values will usually be the same, but you may need to set Reallow as something else.

Link atoms:

If the two layers are not connected through covalent bonds, this issue does not apply. If they are, then link atoms must be generated to fill out the valencies in Modelhigh (otherwise that level would have dangling bonds). If the calculation is set up in Gaussview, these link atoms are automatically generated as H atoms at an automatically-calculated distance from the Modelhigh atom to which it is bound. These distances can be changed by adding scale factors to the input (see above ref. and/or the manual). Therefore, if you treated a tert-butyl group as a MM low layer, it would be treated in the high layer as a group with the sterics of tert-butyl but the electronics of H (so pick your boundary lines carefully).

ONIOM input file:

Inspect the .com file. You will notice that ONIOM inputs automatically get put in Cartesian coordinates (this is set up to work with MM fragments since they must be in Cartesian coordinates). In each atom specification line, the first non-blank character is a 0. If you would like to freeze that coordinate, replace the 0 with -1.
Note: If you freeze atoms in MM low layer and you are not using electronic embedding, it is crucial to specify Opt=Quadmacro. This makes the redundant internal coordinates in the high layer communicate properly with the low-layer Cartesian coordinates.

Following the atom coordinates, each line has the character L or H (or M) which designates it as belonging to the low or high (or middle) layers. If link atoms were added, they are listed on the line of the atoms which they replace in the high layer, followed by the atom numbers to which they are bound (H atoms attached to atom 1 and 10 in the example below). A custom scale can be added after that.

Below is the relevant section of an input file for an ONIOM geometry optimization of dimethylperoxide, with the oxygens treated with the B3LYP/cc-pVQZ model chemistry and the methyl groups treated with the HF/3-21G model chemistry:

# opt oniom(b3lyp/cc-pvqz:hf/3-21g)

dimethylperoxide 

0 1 0 1 0 1
O               0   -0.987457    2.258703   -0.004221 H
C               0   -1.716673    3.483521    0.109640 L      H      1
H               0   -1.607312    4.048882   -0.792196 L
H               0   -2.752107    3.268347    0.272360 L
H               0   -1.336238    4.049806    0.933952 L
C               0    1.513799    1.033520   -0.000000 L      H     10
H               0    1.870465    1.537970    0.873624 L
H               0    1.870465    0.024714    0.000054 L
H               0    1.870466    1.537875   -0.873679 L
O               0    0.307454    2.513737    0.019488 H

Note: if MM were used for the low level, the following Gaussview-generated bond order information would be necessary (along with Geom=connectivity in the Route section):

1 2 1.0 10 1.0
2 3 1.0 4 1.0 5 1.0
3
4
5
6 7 1.0 8 1.0 9 1.0 10 1.0
7
8
9
10



G3 Energy Calculations Using Previously-Calculation Structures:

There are four types of G3 calculations available in Gaussian: G3, G3(MP2), G3B3 and G3(MP2)B3. Those jobs containing MP2 in the name use MP2 for the parts of the calculation normally employing MP4, which saves a lot of time at the expense of some of accuracy. By default, all of these jobs do a geometry optimization and frequency calculation before the compound energy calculations begin. Those jobs containing B3 in the name use B3LYP for these parts. If you want to obtain a G3 energy of a transition state, or a minimum calculated using any basis set other than 6-31G(d), you will need to start the G3 job from a checkpoint file from a frequency calculation containing a geometry and force constant matrix (Hessian) computed previously at your preferred level of theory.

Run your optimization (or TS search) and frequency calculation (opt and freq keywords). When it finishes, make a copy of the checkpoint file, naming it what you will name your G3 calculation. For example, if you name your new checkpoint file g3test.chk, you would then name the new input file g3test.com.

To make the input file, you could either read the new checkpoint file into Gaussview and set up the file there, or create a file using a text editor containing the following:

%chk=/home/username/dirname/filename.chk
%mem=1GB
%nproc=1
# g3(startmp2) geom=allcheck
{blank line}

The order of lines is important here. The %check does not need to contain the full path as long as you submit the job from the same directory containing the checkpoint file. There are only two necessary keywords in the Route section. The keyword g3(startmp2) specifies a G3 calculation with NO geometry optimization (you could use g3mp2, g3b3 or g3mp2b3 here instead of g3).
Note: if you used DFT for your optimization, you still need to specify g3b3 or g3mp2b3 for the G3 single point calculation, or it will still do an MP2 geometry optimization.

No basis set is specified because the compound energy calculation uses a variety of pre-defined basis sets. The other command is geom=allcheck, which calls for all molecule information to be read from the checkpoint file. No molecule information (charge, atomic coordinates, etc) is necessary in the input file. There should be a blank line after the route section.


NBO (Natural Bond Orbital) Analysis:

NBO is a separate program but it is built into our copy of Gaussian. It takes a QM wavefunction and converts it into localized orbitals analogous to valence bond theory. Many features are available but this section only concerns the generation of basic output and visualization of the NBOs. All of the text-based output of an NBO analysis is printed to filename.out. For more information on NBOs see the NBO3 manual (www.ccl.net/cca/software/MS-WIN95-NT/mopac6/nbo.htm) or the NBO website (nbo6.chem.wisc.edu/).
Note: because our copy of Gaussian has NBO version 3.1 built in, it impossible for us to upgrade NBO to a more current version at this time. Thus, do not use the NBO5 documentation as a reference for calculation setup. However we do have NBOView version 5 which can accept plotting files generated in NBO 3.1. So you can use NBOView or the steps below for viewing NBO data.

This how-to assumes you have already done a QM calculation on your molecule. It is a good idea to do the NBO analysis in a separate step because the NBOs will replace any MOs present in the checkpoint file (you probably also want MOs).

Start by copying the checkpoint file from a previous calculation on the molecule of interest (must include an SCF calculation). Read it into Gaussview, and in the Gaussian Setup window select Job Type: Energy, NBO Type: Full NBO, Checkpoint Save: Save NBOs. In the Additional Keywords field, type guess(read,only). If you are setting up the calculation in a text editor, your Route section should now look like this:

# Hamiltonian/Basis pop=(nbo,savenbo) guess=(read,only)

where Hamiltonian and Basis are whatever model chemistry you are using. The guess options make the program read the wavefunction from the checkpoint file and skip the SCF step in the new calculation, making the NBO calculation run very quickly. The savenbo argument replaces the MOs in the checkpoint file with NBOs, so if you are just interested in the text output you do not need that.

When the job completes, read the checkpoint file into Gaussview and visualize the NBOs in the Molecular Orbital Editor. Note: the orbital numbers in the NBO analysis report (in filename.out) do not correspond to those in Gaussview. This is because Gaussview orders them by energy, whereas NBO organizes them by something relating to atom number.


Obtaining Average Quality Graphics Files:

You have two options for obtaining pictures from Gaussview.

1) You can export the molecular model in the current window by saving the file in mol2 or pdb format. This file can be imported for viewing in various visualization programs such as Mercury, Molden, Molekel (see below), VMD or even ORTEP.


2) To obtain a picture directly from Gaussview, you must use an external image capture program (the "Save Image" feature in Gaussview does not function with the new Linux operating system). First you can set your desired background color by selecting Preferences from the File menu in Gaussview, and then selecting Colors from the left of the menu.

Next, open a new command line terminal. Make sure the Gaussview window containing your molecule is not partially obscured by any other windows and your molecule is as large as desired (larger picture = better image quality). In your new terminal, type the command import filename.extension, where filename is whatever you call it and extension is the type of file you wish to create (jpg, tif, bmp, etc.). After executing this command, left click anywhere in the Gaussview window containing your molecule. A beep will be heard and your image file will be saved to the directory from which the command was executed.

To view the picture, type display filename.extension and the file will open in ImageMagick. Clicking anywhere in the window will bring up the control panel which can be used to edit and resave the image.


High-Quality MO (or NBO) Pictures:

If photo-like pictures are desired, it is best to generate cube files in Gaussview and import them into a different program. We used to recommend Molekel (see below) but now we recommend VMD. Click on the links for program description and instructions for making MO pictures with VMD.

For Molekel users, under the Results menu in Gaussview, select Surfaces. This brings up the Surfaces and Cubes panel. Select Cube Actions: New Cube and select what surfaces you want to generate (use Grid = Fine for smooth surfaces). Click Ok, and wait for text to appear in the Cubes Available field. When the surface you want is in the field, select Cube Actions: Save Cube (explicitly type the .cube extension in the filename).

Open Molekel (type molekel on a command line) and open filename.cube. In the white space to the left of the screen, click on filename.cube (a white box appears around the molecule). Under the Surfaces menu, select Grid Data. You must now enter the Value (isovalue), Color, and Rendering Style for your surface. The optimum value for Value varies greatly, so start with the default and experiment from there.

Currently, only one phase of an MO can be generated. A workaround is to import the same structure twice and calculate the same surface for each but with opposite signs for Value. When doing this, make sure you are in Interact with Camera mode (default). If you use Interact with Molecule, you will no longer have superimposed entries.

Under the Display menu, select Molecule to change the molecular representation, then save the picture in your preferred format.

Back to the Graphics Facility Home Page