A typical set of GAMESS calculations

In this post I show the sequence of GAMESS calculations necessary to compute the free energy of the water dimer molecule at the B3LYP/6-31G(d)//PM3 level of theory. This notation means that the geometry and the free energy contribution is obtained with PM3, while the electronic energy is computed using B3LYP/6-31G(d) and the PM3 geometry. (By default the free energy is computed at 298.15 K and 1 atmosphere pressure.)

The steps are as follows:

1. Build the molecule

2. Optimize the geometry using PM3

3. Compute the frequencies for the optimized geometry to
a. verify that you found a minimum
b. compute the sum of the translational, rotational, and vibrational free energies at the PM3 level.

4. Compute the B3LYP/6-31G(d) energy for the PM3 optimized geometry

Some of the keywords I used have been described in a previous post, and I actually copy the keywords from this post to save typing.

The electronic energy for the water dimer is -152.7533 atomic units (au) while the sum of the translational, rotational, and vibrational free energies is 13.820 kcal/mol. If you repeat these calculations for the water molecule you get -76.3715 au and 2.271 kcal/mol (1 au = 627.51 kcal/mol).

So, the change in electronic energy on going from two water molecules to the water dimer is

ΔE = (-152.7533 - 2(-76.3715))*627.51 = -6.4 kcal/mol

The corresponding free energy change is

ΔG(298K) = -6.4 + 13.820 - 2(2.271) = -6.4 + 9.3 = 2.9 kcal/mol

Here I make the usual assumption that the electronic free energy is the electronic ground state energy.

Notice that the free energy change is positive meaning the water dimer is not predicted to form at this temperature and pressure (1 bar) under equilibrium conditions (the relative probability of observing the water dimer can be computed with this equation). This is because of the entropy of loss on going from 2 particles to 1.

Some GAMESS input basics

Update 2 August 2009: I have changed the discussion regarding the SG1 keyword. New text is in italics.

I have planned some posts in which some details of the GAMESS input will become important, so this post is about the basic anatomy of a GAMESS input file, and some non-default settings that I use all the time.

Here is an input file for a single point energy calculation at the B3LYP/6-31G(d) level of theory for the H₂ molecule:

$contrl runtyp=energy dfttyp=b3lyp $end
$basis gbasis=n31 ngauss=6 ndfunc=1 $end
$data
Title line: you can write anything here
c1
H 1.0 0.0 0.0 0.0
H 1.0 0.753 0.0 0.0
$end


• The order of the groups doesn’t matter, so you can put the $basis group below the $data group.

• You can have several groups of the same kind (e.g. two $contrl groups). If you have the same keyword in both groups, it is the value in the last group that counts. For example,

$contrl runtyp=energy dfttyp=b3lyp $end
$contrl runtyp=optimize $end

is perfectly legal, and will result in an geometry optimization.

• The order of the keywords within a group doesn’t matter.

$contrl runtyp=energy dfttyp=b3lyp $end

is the same as

$contrl dfttyp=b3lyp runtyp=energy $end

• Note that there is a space in front of any line starting with a $ sign. This is important and if you leave it out GAMESS will not read that line. Note you can use this to “uncomment” a line by adding, for example, a ! For example, this

$contrl runtyp=energy dfttyp=b3lyp $end
!$contrl runtyp=optimize $end

will result in a single point energy calculation

• You can use capital letters if you want.

• For parameters that can be true for false, you can use either “.true.” or “.t.” and “.false.” or “.f.” (there's an example below).

• A line can only be 80 characters long. Anything after that is ignored. It is perfectly legal to split a line, for example:

$contrl runtyp=energy
dfttyp=b3lyp $end

Some useful options

Almost all of the parameters have default (but there is no default basis set and no defaults in $data). For example, the default runtyp is energy, so

$contrl dfttyp=b3lyp $end

will result in a single point energy evaluation. If $contrl is missing entirely, GAMESS will compute a single point energy at the RHF level.

There is a few defaults that I almost always change

Single point energy calculations (not AM1 or PM3)

$scf dirscf=.t. $end

This tells GAMESS to recompute the integrals at each SCF step rather than reading them from hard disk. Most computers can recompute a number faster than they can find it on the hard disk.

DFT single point energy calculations

$dft sg1=.t. $end

Computing the DFT energy requires a numerical integration over a grid. The default grid is very fine (has a lot of points). Standard Grid 1 (sg1) is adequate for most purposes and much faster. Note that this leads to inaccurate gradients and should not be used for geometry optimizations.

DFT, AM1, and PM3 frequency calculations

$force nvib=2 $end

DFT and semiempirical Hessians (second energy derivatives) are computed numerically, by displacing each atom in the x, y, and z direction and computing the first derivative. This is usually not accurate enough. nvib=2 leads to a displacement in both the positive and negative x, y, and z direction resulting in more accurate frequencies.

Geometry optimizations

$statpt opttol=0.0005 nstep=50 $end
$contrl nzvar=1 $end
$zmat dlc=.t. auto=.t. $end

I think the default criterion for geometry convergence (0.0001) is too strict, and the default number of steps (2o) is too small. So I usually use 0.0005 and 50, respectively. Optimizing in internal coordinates instead of Cartesian coordinates (default) usually speeds up convergence, because the former are "more quadratic". nzvar=1 tells GAMESS to use internal coordinates and dlc and auto tells GAMESS to automatically generate delocalized internal coordinates.

Unfortunately, the DLC implementation in GAMESS only looks for covalent bonds, not H-bonds and such. So you have to specify these if you have specify non-covalently bonded molecules in the input file. For example, in the case of water dimer (see figure for numbering), you do it like this

$zmat nonvdw(1)=1,4 $end

Note that you don't have to specify the H-bond in HO-CH2-CH2-OH, since the OH groups are connected by covalent bonds.

So, an input file for a geometry optimization of the water dimer would look like what I show below. That might seem a lot to type, but the trick is of course to type is once, store it somewhere handy, and copy it into new files.

$contrl runtyp=optimize dfttyp=b3lyp $end
$basis gbasis=n31 ngauss=6 ndfunc=1 $end
$scf dirscf =.t. $end
$statpt opttol=0.0005 nstep=50 $end
$contrl nzvar=1 $end
$zmat dlc=.t. auto=.t. nonvdw(1)=1,4 $end
$DATA
Title
C1
O 8.0 -4.42170 -0.11056 -2.11879
H 1.0 -3.53334 -0.03874 -1.70405
H 1.0 -4.75257 0.79398 -2.01675
O 8.0 -2.06207 0.35222 -0.78539
H 1.0 -1.13417 0.41136 -1.06880
H 1.0 -1.98676 0.07447 0.14406
$END

Installing GAMESS on a Mac

NUchem (could be a pseudonym) suggested a screencast about installing GAMESS. In my case this will have to be installation on a Mac, and the installation on a Windows or Linux computer will be quite different.

If others want to volunteer similar screencasts of Windows or Linux do leave a comment on the blog. I'd be happy to post them or provide a link.

I use GAMESSQ to run GAMESS. It is also possible to run GAMESS from the command line, but I find GAMESSQ much more convenient, so this is what I suggest you use on a mac.

Building a complicated molecule: 2D to 3D

In a previous post I showed how to build a complicated molecule with Avogadro. Here I show how to build it in 2D first and then convert it to 3D.

I use PubChem Editor to build a 2D structure and insert it into Avogadro using the insert SMILES option. JChemPaint is another free option for drawing the 2D structure.

Avogadro does not process the chirality info (the @ signs) in the SMILES string correctly, so some more editing is needed. Hopefully this is fixed in future versions of Avogadro (I use 0.9.5).

There is some discussion over at Rich Apodaca's blog about whether building in 2D or 3D is easier. Try for yourself and make up your own mind.

Update: the program Marvin Sketch can also be used for 2D to 3D conversion, and the chirality info is interpreted correctly.

Building a complicated molecule and finding the lowest energy conformer

Henry Rzepa made an interesting post on rationalizing the difference in half lives of amide hydrolysis in closely related molecules using molecular mechanics. Why not try it for yourself using Avogadro?

The screencast shows you how to build the molecules (using the autoopt tool) and features the automated conformational search option in Avogadro. Enjoy!

The mysterious GAMESSQ program

A recent question from Jacob Lerche mentioned a queueing program for GAMESS called GAMESSQ. Never heard of it. A google search only revealed that the latest MacMolPlt version has a "hook to it". Nothing on the GAMESS page. Strange ...

A few days later the another google search listed Jacob's question, that I posted on the blog comments, as one of the top hits, so now Molecular Modeling Basics looked like an authority on the subject and I decided to look into it.

The program is indeed bundled into the latest GAMESS download (which you can initiate here). It came with an brief manual in html format, so I put it on my web server here. The "about GAMESSQ" feature on the Mac tells me that the program is open source and written by one Jason Ekstrand and "commissioned" by MacMolPlt author Brett Bode.

Anyway, it is a very useful program, and I made a screencast about how to use it. What I don't mention in the screencast (which is already quite long) is (1) that the Jobs>open in wxMacMolPlt doesn't seem to work and (2) the program allows you to submit several jobs at once, which will then run one-after-another. Which of course is the whole point of a queueing program.

10 August, 2009. Update: GAMESSQ has now been officially released.

The AutoOpt tool in Avogadro

I am have been meaning to do a screencast of the Auto Optimization feature in Avogadro for a while. Fortunately Geoff Hutchison beat me to it, as you can see above, and did a much better job than I could have. To prove that point I made the screencast below. It shows how the autoopt tool can be used to illustrate intermolecular interactions - in this case hydrogen bonding between water molecules. Note how the H atoms on one water follow the O atom on the other when you move it (electrostatic attraction), how it's impossible to get the O atoms next to each other (electrostatic repulsion), and how one molecule gets out of the way when you try to push the other too close (steric repulsion).

The ammonia VSEPR

The H-N-H angle in NH₃ is 107^o while the H-C-H angle in CH₄ is 109.5^o. The simplest explanation for both is the Valence Shell Electron Pair Repulsion (VSEPR) theory. Methane adopts a tetrahedral geometry to minimize repulsion between the four electron pairs in the C-H bonds. Similarly for ammonia, but it has a lone pair, which is fatter near the nucleus, so the lone pair-bond repulsion is slightly larger than the bond-bond repulsion. This results in a H-N-H angle that is slightly smaller than 109.5^o.

Here is a Jmol application I wrote to illustrate just that. I use GAMESS and B3LYP/6-31G(d) to compute the geometry (H-N-H angle = 106^o) and molecular orbitals. Molecular orbitals are not unique and most quantum programs do not produce VSEPR orbitals by default. Instead you have to tell the programs to localize the orbitals after they have been computed. There are several criteria for what localized orbitals are, and one of them is to minimize electron pair repulsion as VSEPR theory suggests. In GAMESS (the only program I know of which this option) such localized molecular orbitals are called Ruedenberg orbitals after Klaus Ruedenberg who first devised an algorithm to compute them with modern theoretical methods.

Jmol allows you to load several files at once and display orbitals for each simultaneously. So I simply computed orbitals for three different ammonia orientations (get the output files here, here, and here). Then it was simply a matter of writing a small script to go from one geometry to the other. It's very simple and the entire script can be found here:

load files "nh3b.log" "nh3c.log" "nh3d.log"
background white
rotate x -90
frame 1.1;isosurface mo26 mo 26 fill translucent 0.4
delay 4
frame 1.1;isosurface mo24a mo 24 fill translucent 0.9;delay 0.2
frame 1.1;isosurface mo24a mo 24 fill translucent 0.8;delay 0.2
frame 1.1;isosurface mo24a mo 24 fill translucent 0.7;delay 0.2
frame 1.1;isosurface mo24a mo 24 fill translucent 0.6;delay 0.2
frame 1.1;isosurface mo24a mo 24 fill translucent 0.5
delay 3
frame 1.1;isosurface mo24a mo 24 fill translucent 1.0
frame 3.1;isosurface mo24d mo 24 fill translucent 0.5;
select 3.1; spacefill off; wireframe off;
select 2.1; spacefill off; wireframe off;
frame all
delay 1
frame 3.1;isosurface mo24d mo 24 fill translucent 1.0;
frame 2.1;isosurface mo24 mo 24 fill translucent 0.5
frame all
spin on

You know, for kids!

Oxygen from Christopher Hendryx on Vimeo.

Much of this blog is about making chemistry more accessible through animation. While animation helps greatly in teaching at the college level, it is absolutely crucial when trying to connect with the general public. Creating such animations is also much more difficult since it has to be entertaining in addition to being informative.

So hats off to Christopher Hendrix for doing just that. This movie is a great service to chemistry. Now it's up to us to help bringing it to the attention of as many people as possible. Take a minute to think about a web page, family member, facebook friend, congressman who could benefit from a link to this movie ... good! Now, why not send them an email?