How to run GAMESS-UK calculations
This document will show you how to set up and run a GAMESS-UK calculation from scratch. It will cover building the initial structure, generating the GAMESS-UK input file, running a calculation, performing an occupation analysis on the molecular orbitals and producing images of the HOMO and LUMO. Additional information on most of the concepts introduced in this documant can be found in the GAMESS-UK user manual and the EaStCHEM Utilities user manual.
In this tutorial we are going to perform a geometry optimization on catechol and then examine the nature of the HOMO and the LUMO through both numerical and graphical analyses.
Building the Initial Structure
There are many different ways of generating a starting structure for a molecular electronic structure calculation. Here, I am going to use ChemSketch to generate a MOL file from which we can start our geometry optimization. I simply draw the molecule; click '3D optimization' and then '3D viewer'. This gives me a 3D representation of catechol with reasonably sensible values for the bond lengths and angles. I then save this as a MDL MOL file (catechol.mol).
Generating the GAMESS-UK Input
We can generate a starting point for our GAMESS-UK input file directly from out MOL file by using the mol2inp utility on the clusters. To start with, enter the following command:
mol2inp -o gam catechol.mol catechol.in
The -o gam option tells it to produce input for GAMESS-UK. This will generate a catechol.in file that looks like
title Created by mol2inp geometry angstrom 0.8855 -0.4840 0.6163 6.0 C -0.3159 -0.7690 1.3869 6.0 C 0.7904 0.2583 -0.6320 6.0 C -1.6124 -0.3116 0.9089 6.0 C -0.5060 0.7156 -1.1102 6.0 C -1.7075 0.4305 -0.3396 6.0 C 1.9697 0.5377 -1.3876 8.0 O 2.1581 -0.9328 1.0846 8.0 O -0.2397 -1.3639 2.3877 1.0 H -2.5755 -0.5399 1.5266 1.0 H -0.5821 1.3106 -2.1109 1.0 H -2.7468 0.7969 -0.7228 1.0 H 1.9007 1.0796 -2.2992 1.0 H 3.0351 -0.7246 0.5216 1.0 H end basis 3-21G runtype scf scftype rhf enter
We will now go through the input file line-by-line explaining the meaning. This will give us enough understanding to be able to alter the file to include the options we want.
The first two lines give our job a name.
title Created by mol2inp
The next block specifies the starting geometry for our calculation.
geometry angstrom 0.8855 -0.4840 0.6163 6.0 C ... ... 3.0351 -0.7246 0.5216 1.0 H end
The line geometry angstrom tells GAMESS-UK that the input geometry is in Cartesian coordinates (rather than a z-matrix) and that the units are Angstroms. Then there is a line for each atom in the molecule with its x, y, and z coordinates, its atomic number and its label. The end line finishes the geometry input block.
The next line,
basis 3-21G
specifies the basis set to use for the calculation. This specification can be as simple as one line (as here) or much more complex, specifying a different basis set for each atom and even specifying a completely custom basis set, GAMESS-UK user manual has more details.
Next we specify the runtype.
runtype scf
This is the type of calculation that we want to perform. There are a number of different types available but the most common are scf (energy at a fixed geometry), optxyz (geometry optimization), analyse (wavefunction analysis) and hessian (vibrational analysis.)
The line
scftype rhf
specifies that we want to use the RHF (Restricted Hartree-Fock) method for the calculation. Other valid options would be uhf (for unrestricted Hartree-Fock) or mp2 (for Moller-Plesset perturbation theory.)
The final line
enter
tells GAMESS_UK to start the calculation using the options specified above.
As we are going to run a DFT geometry optimization (using the B88 exchange and LYP correlation functionals), with a 6-31G(d) basis set we will have to change the input file. Open the file with your favourite text editor (either use one on the cluster or download it to your PC and use Notepad) and change it to look like (catechol.in)
file ed3 catechol.ed3 title Created by mol2inp geometry angstrom 0.8855 -0.4840 0.6163 6.0 C -0.3159 -0.7690 1.3869 6.0 C 0.7904 0.2583 -0.6320 6.0 C -1.6124 -0.3116 0.9089 6.0 C -0.5060 0.7156 -1.1102 6.0 C -1.7075 0.4305 -0.3396 6.0 C 1.9697 0.5377 -1.3876 8.0 O 2.1581 -0.9328 1.0846 8.0 O -0.2397 -1.3639 2.3877 1.0 H -2.5755 -0.5399 1.5266 1.0 H -0.5821 1.3106 -2.1109 1.0 H -2.7468 0.7969 -0.7228 1.0 H 1.9007 1.0796 -2.2992 1.0 H 3.0351 -0.7246 0.5216 1.0 H end basis sv H 6-31G svp C 6-31G svp O 6-31G end runtype optxyz scftype rhf dft B88_X dft LYP punch coordinates basis vectors occupation eigenvalues scfenergy 1 enter
We have added a few lines
file ed3 catechol.ed3 - This line tells GAMESS-UK to save the dumpfile, this file is needed to restart the calculation if it crashes or we want to run some analysis on the results.
The basis section is now made up of separate basis specifications for each of the atom types. sv H 6-31G specifies a split-valence basis set (6-31G) for hydrogen while svp C 6-31G and svp O 6-31G specify split-valence basis sets with polarization functions for carbon and oxygen. The basis specification block is terminated by end. This basis block defines the 6-31G(d) basis set for the molecule.
runtype optxyz - Tells GAMESS-UK to run a geometry optimization.
dft B88_X - Specifies the B88 exchange functional.
dft LYP - Specifies the LYP correlation functional.
punch coordinates ... - Tells GAMESS-UK to write out details of the orbitals to the punch file. The punch file is used for storing formatted output from GAMESS-UK and is used to calculate orbital occupancies. We will see later that it can also be used to plot the orbitals. The final integer on the line specifies the section of the dumpfile to et the orbitals from. For RHF and DFT calculations, this is always 1.
Now that we have our input we will move on to running the calculation.
Running the Calculation
We are going to use the EaStCHEM Software Submission Package (SSP) to submit our job to the batch queuing system. This package automatically generates our submission script for us.
To submit our job we use the gamesssub command. Type the following on the command line:
gamesssub -np 4 -q parallel-short.q catechol
This should result in output that looks something like this
Creating GAMESS-UK job:
Local job on alchemy.epcc.ed.ac.uk
Using 4 distributed memory processors (4 each on 1 node[s])
Charging to account: rcf
Job File: /users/aturner/test/gamess/catechol/catechol.bash
Submitting local job...
Your job 10801 ("catechol_gam") has been submitted
job-ID prior name user state submit/start at queue ------------------------------------------------------------------------------------------------
10793 0.51484 triplet__g mpalmer r 10/02/2007 11:39:12 fat-medium.q@comp02.epcc.ed.ac
10785 0.55211 PbF2_flu3 dmarroc r 10/02/2007 03:29:28 parallel-medium.q@comp09.epcc.
10641 0.56113 PbF2_cot4 dmarroc r 10/01/2007 14:58:01 parallel-medium.q@comp11.epcc.
that is the job submitted. You now need to wait for it complete before continuing.
Checking the output
Once your calculation is finished, the working directory should contain the following files
[filelist]
- e, o, pe, po are produced by the queuing system and are not needed if the calculation was successful. They can be useful to diagnose errors if things go wrong.
- catechol.out contains the main output from the GAMESS-UK calculation.
- catechol.ed3 is the dumpfile that we specified at the top of our input.
- catechol.pun is the punch file containing all the information we asked for in a well defined format.
To check that the calculation completed successfully we need to look at the catechol.out file. We will do this using less (a command that prints the contents of a file to the screen). Type
less catechol.out
This should display the top of the catechol.out file. While using less you can use the following keys: g - go to top of file, G - go to bottom of file, space - down a page, u - up a page, return - down a line, q - quit less and return to command line. We are going to use another function of less, th search function. Type the following (while using less)
/optimisation converged
The slash tells less that we want to search for the specified string in the file. If you calculation has worked you should see the following
[output]
Now we have checked that everything is OK we can have a closer look at the output file to identify useful information. One of the most basic pieces of information we can get from the calculation is the total molecular energy. This property allows us, for example, to identify the most stable conformer of a molecule or compute binding energies.
Calculating the Orbital Positions
It is often very useful to know how much of a particular MO is located on a particular atom (or basis function.) Using the MO coefficients from the punch file written by the above calculation we can calculate this (it is computed from the square of the orbital coefficients - the orbital coefficients actually specify how much of a given basis function is in that MO.) Performing the analysis is simple using the getocc utility. Type the following command
getocc -c gam catechol.pun catechol.mo
The -c gam option tells the utility that the input is a GAMESS-UK punch file (catechol.pun). THe output from the analysis will be in catechol.mo. Now we need to take a look at the catechol.mo file using the less command:
bash> less catechol.mo
**********************************
* getocc *
* MO Occupancy analysis *
* A. R. Turner, EaStCHEM RCF 2007 *
***********************************
Molecular Orbitals
Input file: catechol.pun
===============================================================================
MO 1, Energy = -18.83790 a.u., Occupancy = 2.00
===============================================================================
By basis function: By atom:
# Atom Func c c^2 % # Atom %
--------------------------------------------- ---------------------
1 c 1 s -0.00002 0.00000 0.00 1 c 0.00
2 c 1 s -0.00014 0.00000 0.00 2 c 0.00
3 c 1 p(x) -0.00006 0.00000 0.00 3 c 0.00
4 c 1 p(y) 0.00001 0.00000 0.00 4 c 0.00
5 c 1 p(z) 0.00000 0.00000 0.00 5 c 0.00
catechol.mo lines 1-23/18754 0%
We are really interested in the HOMO and LUMO but the file contains the analysis for all the MOs. How do we identify the HOMO and LUMO? The answer is that we need to know how many electrons are in the system and divide it by two (as each orbital is doubly occupied in a restricted DFT calculation.) NB. If we had performed an unrestricted calculation, used an ECP basis set or a calculation with fractional orbital occupancies then this process would be more complicated. We can either add up the electrons in the system manually or check the GAMESS-UK output. I am going to use grep (a command that searches text files for the specified string and prints lines where it occurs):
grep "number of electrons" catechol.out number of electrons 58
We have 58 electrons so the HOMO is MO number 29 and the LUMO is MO number 30. If we go back to the catechol.mo file now we can search for the MOs we are interested in (as we did to check if the optimisation was converged) by typing / 29, (space-29-comma):
bash> less catechol.mo
**********************************
* getocc *
* MO Occupancy analysis *
* A. R. Turner, EaStCHEM RCF 2007 *
***********************************
Molecular Orbitals
Input file: catechol.pun
===============================================================================
MO 1, Energy = -18.83790 a.u., Occupancy = 2.00
===============================================================================
By basis function: By atom:
# Atom Func c c^2 % # Atom %
--------------------------------------------- ---------------------
1 c 1 s -0.00002 0.00000 0.00 1 c 0.00
2 c 1 s -0.00014 0.00000 0.00 2 c 0.00
3 c 1 p(x) -0.00006 0.00000 0.00 3 c 0.00
4 c 1 p(y) 0.00001 0.00000 0.00 4 c 0.00
5 c 1 p(z) 0.00000 0.00000 0.00 5 c 0.00
/ 29,
this will skip down the file to display
MO 29, Energy = -0.16734 a.u., Occupancy = 2.00
===============================================================================
By basis function: By atom:
# Atom Func c c^2 % # Atom %
--------------------------------------------- ---------------------
1 c 1 s 0.00000 0.00000 0.00 1 c 15.18
2 c 1 s -0.00000 0.00000 0.00 2 c 3.79
3 c 1 p(x) -0.00002 0.00000 0.00 3 c 15.73
4 c 1 p(y) -0.00001 0.00000 0.00 4 c 9.33
5 c 1 p(z) -0.27751 0.07701 11.09 5 c 0.48
6 c 1 s 0.00002 0.00000 0.00 6 c 15.49
7 c 1 p(x) -0.00002 0.00000 0.00 7 o 15.96
8 c 1 p(y) -0.00002 0.00000 0.00 8 o 24.05
9 c 1 p(z) -0.16710 0.02792 4.02 9 h 0.00
10 c 1 d(xx) 0.00000 0.00000 0.00 10 h 0.00
11 c 1 d(yy) -0.00000 0.00000 0.00 11 h 0.00
12 c 1 d(zz) 0.00000 0.00000 0.00 12 h 0.00
13 c 1 d(xy) 0.00000 0.00000 0.00 13 h 0.00
14 c 1 d(xz) 0.01533 0.00024 0.03 14 h 0.00
15 c 1 d(yz) -0.01568 0.00025 0.04 ---------------------
16 c 2 s -0.00000 0.00000 0.00 Total 100.00
17 c 2 s 0.00001 0.00000 0.00
catechol.mo lines 3989-4011/18754 21%
We can see (on the far right) that the HOMO is 24.05% on O8 and 15.96% on O7. If you repeat this for the LUMO (MO 45) you should see that the LUMO is 31.17% on C5.
This is all well and good but in practice we often want to see a graphical representation of the orbital. In order to do that we must restart the calculation and tell GAMESS-UK to compute 3D representations of the MOs that we are interested in.
Restarting the Calculation
We are going to create a new GAMESS-UK input file that will read the results of the previous optimization from the dump file (catechol.ed3) and use them to produce 3D representations of the HOMO and LUMO.
The input file (catechol_plot.in) for this restart and analysis job looks like this
file ed3 catechol.ed3
restart
title
Catechol graphical analysis
geometry angstrom
0.8855 -0.4840 0.6163 6.0 C
-0.3159 -0.7690 1.3869 6.0 C
0.7904 0.2583 -0.6320 6.0 C
-1.6124 -0.3116 0.9089 6.0 C
-0.5060 0.7156 -1.1102 6.0 C
-1.7075 0.4305 -0.3396 6.0 C
1.9697 0.5377 -1.3876 8.0 O
2.1581 -0.9328 1.0846 8.0 O
-0.2397 -1.3639 2.3877 1.0 H
-2.5755 -0.5399 1.5266 1.0 H
-0.5821 1.3106 -2.1109 1.0 H
-2.7468 0.7969 -0.7228 1.0 H
1.9007 1.0796 -2.2992 1.0 H
3.0351 -0.7246 0.5216 1.0 H
end
basis
sv H 6-31G
svp C 6-31G
svp O 6-31G
end
runtype analyse
punch coordinates basis grid 160 161
graphics
gdef
title
3D Grid
type 3d
points 60
size 15.0
calc
title
MO 29 (HOMO)
type mo 29
section 160
calc
title
MO 30 (LUMO)
type mo 30
section 161
enter
There are a lot of changes from our original calculation input so I will take you through them one-by-one.
The addition of the restart directive before the title tells GAMESS-UK to use the details from the ed3 file.
- We still have to include a geometry section even though the geometry is taken from the ed3 file. Do not know why!
It is now runtype analyse as we just want to print out the orbitals.
The punch directive includes the options for printing the orbitals to the formatted punch file.
We use grid 160 161 to print the orbital plots stored on sections 160 and 161 of the ed3 file (see below.)
The section starting with graphics defines the orbitals we want to print:
First we must define a 3D grid to plot our orbitals on to. This is what the gdef section does. The only interesting option is the size 15.0. This defines the grid as being a 15 angstrom grid centred on the molecule. The size should be big enough to encompass your whole molecule.
The calc sections define the orbitals to plot. All you need is a title, type mo 29 specifies we want MO #29 and section 160 tells GAMESS-UK to write the plotted orbital to section 160 of the ed3 file.
Repeat the calc section for the LUMO but make sure it writes to a different section.
OK, now we can run this calculation. Of course, we have to make sure that our old ed3 file is accessible to the new calculation. This is a non-standard running scenario so we will have to alter our job submission file before submitting the job. Fortunately this is relatively easy as gamesssub allows us to create a template to work from. To create a template rather than actually run the job we just add the -norun option.
gamesssub -norun -q serial-medium.q catechol_plot
This will have created a template job submission script called catechol_plot.bash (for more on submission scripts see How To Write SGE Submission Scripts.) As it is just to plot the orbitals and not do any serious calculatio we can just use a serial queue rather than the parallel queue we used for the actual calculation. Open this file using a text editor and add the extra line as shown below
#!/bin/bash # #$ -cwd -V #$ -q serial-medium.q #$ -N catechol_plot_gam #$ -A rcf #$ -l medprios=true -l h_rt=03:00:00 cwd=`pwd` cd $scratch cp $cwd/catechol.ed3 $scratch/ # <== This is the line that has been added /usr/local/GAMESS-UK-7.0/bin/gamess < $cwd/catechol_plot.in > $cwd/catechol_plot .out if [ -e ftn058 ]; then mv ftn058 $cwd/catechol_plot.pun fi if [ -e catechol_plot.ed3 ]; then mv catechol_plot.ed3 $cwd/catechol_plot.ed3 fi
This additional line copies the ed3 file from our previous calculation into the directory that our plotting calculation will run in.
Now we can submit the job to the queue system
> qsub catechol_plot.bash
Once the calculation has finished your directory should contain the following files
[ls output]
Next we will see how you can view the orbitals in 3D and produce pictures for publication.
Viewing the Orbitals
Although GAMESS-UK provides its own visualiser to view orbitals and such it is hampered by being rather unwieldy and has an inability to produce high quality pictures for publication. To overcome these limitations we are going to convert our MO plots to the Gaussian Cube file format. This will allow us to view the orbitals using more convenient programs and also produce high quality images.
To do this we will use the pun2cub utility on the clusters. Firstly, we will check that we have the orbitals in the punch file, type
> pun2cub -list catechol_plot.pun
and you should see the following output
Input File = catechol_plot.pun
block title
===== =====
1 MO 29 (HOMO)
2 MO 30 (LUMO)
This tells you that there are two plotting blocks in the punch file: block #1 is the HOMO and block #2 is the LUMO. We will no extract the HOMO along with a molecular structure that the orbitals can be plotted on to
> pun2cub -block 1 -type mo -xyz catechol_opt.xyz -cube catechol_homo.cub catechol_plot.pun
Here we are extracting block #1 (the HOMO) it is of type mo (a MO), we will write the structure to catechol_opt.xyz and the cube file to catechol_homo.cub. We will now repeat this for the LUMO (block #2) but without producing the structure file (as we already have it!)
> pun2cub -block 2 -type mo -cube catechol_lumo.cub catechol_plot.pun
Now we need to transfer the two .cub files and the catechol_opt.xyz file to our PC so we can view them.
We are going to use a free program called Argus Lab to view our MOs You can download it from. To start open Argus Lab, click Open and open the catechol_opt.xyz file. This should give you a window with a view of the optimized catechol structure.
[Picture of window]
To plot the HOMO go to Surfaces>Make Surfaces. In the resulting dialog click Import Grid. Locate the 'catechol_homo.cub file and load it. Now drag the picture of the blue cube to the, change the name to HOMO and click create. The dialog should now look like: [Picture of dialog]
Under Construction