Monday, October 31, 2011

Using OpenBabel to Make Input Files for GAMESS

One thing that people tend to hate to do, is to make repetitive tasks. One such thing could be to prepare 100 different geometry optimizations using DFT[1] and a happy mixture of basis sets[2] and functionals[3] for GAMESS[4]. In this blog post, I'll show you how you can get around manual labor and prepare 10 files with just a little bit of coding. Extending it to a 100 is an easy task after that. Lets start.

An input geometry of ethanol in the very common XYZ file format might look something very similar to this, when made in Avogadro[5]

9
ethanol.xyz
C  0.00004  0.00001  0.00002
C  0.00035 -0.00106  1.52091
H  0.02794  1.04324 -0.38022
H -0.92039 -0.49540 -0.37524
H  0.88619 -0.55034 -0.38151
H -0.89666  0.55109  1.87908
O  1.16332  0.62487  1.98684
H -0.03205 -1.05436  1.87773
H  1.11207  0.59582  2.97727

Using Open Babel[6] we can quickly make a GAMESS input file using the following command

babel -ixyz ethanol.xyz -ogamin ethanol.inp

resulting in

$CONTRL COORD=CART UNITS=ANGS$END
$DATA ethanol.xyz C1 C 6.0 0.0000400000 0.0000100000 0.0000200000 C 6.0 0.0003500000 -0.0010600000 1.5209100000 H 1.0 0.0279400000 1.0432400000 -0.3802200000 H 1.0 -0.9203900000 -0.4954000000 -0.3752400000 H 1.0 0.8861900000 -0.5503400000 -0.3815100000 H 1.0 -0.8966600000 0.5510900000 1.8790800000 O 8.0 1.1633200000 0.6248700000 1.9868400000 H 1.0 -0.0320500000 -1.0543600000 1.8777300000 H 1.0 1.1120700000 0.5958200000 2.9772700000$END

but we can't run this file right away since it is obviously lacking information about the basis set, the DFT-functional and even the runtype (GAMESS assumes runtyp=energy if nothing is specified).

Ideally what we want is to generate a header for the input file and then insert it into the input file, but this header should be generated automatically for us.

Given a set of parameters (functional and basis set), the following python script (header.py) prints a correct header for GAMESS to the screen

import sys

basissets = {'sto3g' :'STO NGAUSS=3',
'pc0' :'PC0',
'pc1' :'PC1',
'apc0' :'APC0',
'apc1' :'APC1'}

functionals = {'blyp' :'BLYP',
'b3lyp': 'B3LYP'}

basis = basissets[sys.argv[1]]
functional = functionals[sys.argv[2]]
print " $SYSTEM MWORDS=125$END"
print " $CONTRL RUNTYP=OPTIMIZE DFTTYP=%s ISPHER=1$END" % (functional)
print " $BASIS GBASIS=%s$END" % (basis)

Notice how we use a dictionary to store the basis set information so it can be easily recovered and printed correctly. The GAMESS basis group can be quite hideous with Pople basis sets[7].

All we require now is a small bash script, which makes loops over the various basis sets and functionals, invoke the python script we just made with the correct arguments and dump that printout information to a file (header.inp) and then use some Open Babel magic (to add a custom header with keywords, use the very secret -xf filename flag when converting to GAMESS input files) to make it work. The following script does exactly that

#!/usr/bin/env bash
for basis in sto3g pc0 pc1 apc0 apc1
do
for functional in blyp b3lyp
do
filename=$basis"_"$functional".inp"
python header.py $basis$functional > header.inp
babel -xf header.inp -ixyz ethanol.xyz -ogamin \$filename
done
done

For information on how to use GAMESS, I recommend Professor Jan Jensens blog (http://molecularmodelingbasics.blogspot.com/) as a great resource and for general usage notes of GAMESS and the Google user group (http://groups.google.com/group/gamess) is very helpful if you have any problems or questions related to running GAMESS.

NB! To discover what options the babel file conversion tool supports for your file format (gamin in this case), simply invoke the following command

babel -H gamin

Happy Conversion!