DFTB3 in GAMESS

I recently did my first DFTB3 calculations in GAMESS.  The process wasn't completely obvious, hence this post.

First I got the parameter files from dftb.org.  You have to fill in and email a registration form to get a login so that might take a few days.  I downloaded the 3ob-3-1 parameter set.

The DFTB3-related input can be found below.  My molecule contains C, N, O, and H atoms.  If there are additional atoms then you need to define additional pairs.  It is OK to specify atom pairs for atoms that are not in the molecule, so the easiest might be to define all relevant pairs once and for all and use it in all input files.

The value for dampex and hubder are found in the README file that is located in the 3ob-3-1 folder I downloaded.  hubder defines parameters for N, C, O, and H - in that order.  It was not clear from the manual what the order should be.  I found the correct order by running the input file without hubder. This run prints out the default hubder values in the order  N, C, O, and H and it turns out that is the order you need to use in the input file.

2016.02.09 Update: Yoshio Nishimoto has told me that the order of the hubder parameters is determined by the order in which they appear in $DATA.  Also using MODIO=31 in $SYSTEM will speed up the calculation.

 $basis gbasis=dftb $END
 $dftb ndftb=3 DAMPXH=.t. dampex=4.0
       hubder(1)=-0.1535,-0.1492,-0.1575,-0.1857 $end
 $dftbsk
 C C "/Applications/gamess/3ob-3-1/C-C.skf"
 C H "/Applications/gamess/3ob-3-1/C-H.skf"
 C N "/Applications/gamess/3ob-3-1/C-N.skf"
 C O "/Applications/gamess/3ob-3-1/C-O.skf"
 H C "/Applications/gamess/3ob-3-1/H-C.skf"
 H H "/Applications/gamess/3ob-3-1/H-H.skf"
 H N "/Applications/gamess/3ob-3-1/H-N.skf"
 H O "/Applications/gamess/3ob-3-1/H-O.skf"
 N C "/Applications/gamess/3ob-3-1/N-C.skf"
 N H "/Applications/gamess/3ob-3-1/N-H.skf"
 N N "/Applications/gamess/3ob-3-1/N-N.skf"
 N O "/Applications/gamess/3ob-3-1/N-O.skf"
 O C "/Applications/gamess/3ob-3-1/O-C.skf"
 O H "/Applications/gamess/3ob-3-1/O-H.skf"
 O N "/Applications/gamess/3ob-3-1/O-N.skf"
 O O "/Applications/gamess/3ob-3-1/O-O.skf"
 $end

I would like to thank +Anders Steen Christensen for his help in figuring this out.


This work is licensed under a Creative Commons Attribution 4.0

MOPAC/PCM interface in GAMESS

A new version of GAMESS was released recently.  Among the many new features is +Casper Steinmann's recent work on interfacing the PCM solvation method with the semiempirical methods AM1 and PM3.  So you can now simulate molecules in solution with AM1 and PM3.  

Casper also parallelized AM1 and PM3 (both in gas phase and solution), though molecules have to be pretty big before you see an appreciable speed up, as you can see from the figure below (on small proteins).

The work has been published in the open access journal PLoS ONE. +Jimmy Charnley Kromann is hard at work on implementing PM6.

This work is licensed under a Creative Commons Attribution 3.0 

"Missing" MOPAC parameters

A long, long time ago in a land far, far away I implemented MNDO, AM1, and PM3 in GAMESS.  This was done by taking chunks of code from MOPAC that contained contained the parameters, integral code, and Fock matrix builder.  In doing so, I never noticed that the parameter file contained more parameters than where published in the papers describing the method.  This conundrum came back to haunt us recently as we're trying to implement PM6.

Read the rest of the post over at Proteins and Wave Functions

Installing GAMESS on a linux PC

Probably best viewed in full screen mode

2014.02.18 Update: See this post by +Jimmy Charnley Kromann for more up-to-date instructions

Alchemist requested a post on installing GAMESS on a linux system: after all "I[f] someone can not install a program, how he/she can use it?" Fair enough.

If you have a mac or windows machine, pre-compiled versions of GAMESS is available (see this post on for installing such a version on Macs).  But if you have bought a linux computer, you do have to compile GAMESS yourself.

Three disclaimers:
(1) The post is aimed at people who want to install GAMESS on their personal computer that happens to run linux.
(2) The post is based on installing the March 25, 2010 version of GAMESS and corresponding scripts as they were distributed in mid August, 2010.  If you view this post a few years from now things may have changed.
(3) I don't have a linux PC, so the screencast is made on a shared research computer.  This affects one of the steps as I describe below, and may affect others steps that I don't know about.

The screencast assumes you have already downloaded gamess (gamess-current.tar.Z).  To download GAMESS start here and follow the instructions.  If you have problems watch the first few minutes of the screencast this post.

Here's a summary of the steps in the  screencast:

1. Unpack GAMESS: type "zcat gamess-current.tar.Z | tar -xvf -"
Now you should have a GAMESS directory.  The file gamess/misc/readme.unix has most of the instructions you need.  What follows it basically a condensed version with a few modifications.

2. In the GAMES directory: type "./config".
Follow the instructions and answer the questions.  You need to know whether you have a 32- or 64-bit chip.

3. Type "/sbin/sysctl -w kernel.shmmax=1610612736"
If you are curious about what this does, read section 5 of the file gamess/ddi/readme.ddi.
You need to be logged in as root for this to work.  I don't have root access on the machine I used, so I get an error message in the screencast.  But our system administrator had executed the command previously.

5. Change to the ddi subdirectory (gamess/ddi).  Type "./compddi >& compddi.log" and then "mv ddikick.x .."

6. Go back to the gamess directory and type "./compall >& compall.log"
It will take 30-60 minutes to compile GAMESS.

7. Type  "./lked gamess 00 >& lked.log"

8. Edit the file rungms.  Here are the lines I modify or add:

set SCR=./
set USERSCR=./
set GMSPATH=./

setenv ERICFMT ./ericfmt.dat
setenv MCPPATH ./mcpdata

setenv  MAKEFP $USERSCR/$JOB.efp
setenv   GAMMA $USERSCR/$JOB.gamma
setenv TRAJECT $USERSCR/$JOB.trj
setenv RESTART $USERSCR/$JOB.rst
setenv   INPUT $SCR/$JOB.F05
setenv   PUNCH $USERSCR/$JOB.dat

   if ($os == Linux)   set GMSPATH=./

#  if ($NCPUS == 1) then
#     set NNODES=1
#     set HOSTLIST=(`hostname`)
#  endif
#
   set HOSTLIST=(`hostname`:cpus=$NCPUS)
   set NNODES=1 

#           echo I do not know how to run this node in parallel.
#           exit 20

9. Edit one line in the file runall: #chdir /u1/mike/gamess.
Type "./runall >& runall.log"

10. Go to the tools/checktst subdirectory and type "./checktst"
(Update: if the gamess directory is not in your home directory, you need to change this in the checktst file.  See this comment)

Running GAMESS in parallel
If you computer has more than one core, you may want to run GAMESS in parallel. To run exam01 on 2 cores type "./rungms exam01 00 2 >& exam01.log &"


Other useful programs and getting started with GAMESS
Now that you have GAMESS running you may want to install other useful programs such as Avogadro, GAMESSQ, and MacMolplt.
Here are some related blogposts to wet you appetite:
Using GAMESSQ
Building a complicated molecule with Avogadro

And here are a few blogposts on getting started with GAMESS:
Some GAMESS input basics
A typical set of GAMESS calculations

Acknowledgment: Mike Schmidt and Casper Steinmann helped with this post.

August 30, 2010 update: Alchemist has made the following page on installing GAMESS on Ubuntu 64 bit.  Be sure to check out the comments for additional useful links.

October 20, 2010:  Are you using Ubuntu 9.10 and 10.04, Linux Mint 8 and 9, Knoppix 6.2.1 or Suse 11.2? Be sure to check out this comment by Mott.

Amide hydrolysis, revisited 2

In a previous post I discussed the TSs of acetamide hydrolysis and hydrolysis of a related amide (3).  I have already made a post on how to find the TS for 3, and in this post I summarize the two ways of finding the TS for acetamide hydrolysis, that I describe in Chapter 5 of the book (where you can find many more details).

To start, I use Avogadro to build the geometry shown in Figure 5.32a and use it to construct the input file for an optimization with 4 constrained bond length.  The values for the constraints are taken from a paper.  Figure 5.32b shows the equilibrium geometry obtained with GAMESS, after 17 steps. Based on this geometry GAMESS finds the TS in three steps.

Figure 5.32a. Initial guess geometry for the constrained optimization discussed in Figure 5.33.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

Figure 5.32b. The geometry resulting from the constrained optimization and the normal mode of the imaginary frequency computed for this structure at the PM3 level.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

It will not always be possible to find an article that reports the TS structure for your reaction of interest. So let’s try to find the TS without the information from the article.  I start from the same structure (Figure 5.32a), and I arbitrarily pick C1-O10 as the reaction coordinate constrained to 1.70 Å.

This constrained optimization results (after 23 steps) in the geometry shown in Figure 5.35a.  The subsequent TS search results in the geometry shown in Figure 5.35b, which has no imaginary frequency and is clearly the product.

Figure 5.35a. The geometry resulting from a constrained optimization in which the distance between atoms 1 and 10 (see Figure 5.32a for numbering) was constrained to 1.7 Å, and the normal mode of the imaginary frequency computed for this structure at the PM3 level.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

Figure 5.35b. The geometry resulting from a TS search initiated from the geometry shown in Figure 5.35(a).

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

The fact that the C–O bond is re-formed indicates that it should be stretched more in the initial guess, so I repeat the optimization with the C–O distance constrained to 2.0 Å instead of 1.70 Å. This results (after 39 steps) in the geometry shown in Figure 5.36.  Using this as a starting geometry for the TS search leads to the TS in Figure 4.30a after 17 steps.

Figure 5.36. The geometry resulting from a constrained optimization in which the distance between atoms 1 and 10 (see Figure 5.32a for numbering) was constrained to 2.0 Å, and the normal mode of the imaginary frequency computed for this structure at the PM3 level.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

In the screencast below I try to reproduce the calculations that lead to these figures.  Because I rebuild the structure in Figure 5.32a, I get different starting coordinates, so the energies and number of steps are different than what I discuss in the book for the 4-constraints approach.

In the case of the 1-constraint approach, I actually manage to find the TS using the 1.70 Å constraint.  Thus, the structures in Figures 5.35b and 5.36 do not appear in the screencast.  This just goes to show how finicky TS searchers are to starting geometries.


While editing the screencast I noticed that the PM3 energies of the two TS structures I find are very different (11 kcal/mol).  This is also true for the TSs I found when writing the book (where they are different by 6.6 kcal/mol).

This appears to be a problem that PM3 has with structures where bonds are partially broken or formed.  I could not reproduce this for structures with normal bond lengths such as the reactants and products.   Note that in the book, I use M06/6-31G(d) single point energies to compute barriers, and that the  M06/6-31G(d)//PM3 single point energies of the TSs found with the two different method are only different by 0.03 kcal/mol.

Polarization and intermolecular interaction

Figure 4.16. 0.002 au isodensity surface with superimposed molecular electrostatic potential for (a) methane using the methane density for the methane–water dimer and (b) free methane monomer. The maximum potential value, 0.01 au, is five times smaller than in previous figures to make the increase in negative charge visible. The black spheres in (a) denote the position of the three nuclei of water. The level of theory is M06/6-31+G(2d,p)//M06/6-31G(d).
Click on the picture for an interactive version.
Click here for a pop-up window
From Molecular Modeling Basics CRC Press, 2010

In a previous post I showed how the difference in the strengths of interaction between methane, methane-water, and water dimer is due to their differences in polarity, which can be visualized using molecular electrostatic potentials (MEPs). Methane interacts stronger with water than with methane because of polarization, i.e. rearrangement of the methane electrons due to the polar water molecule that, to a first approximation, induces a dipole moment in methane.

The polarization of the methane density is not visible in the MEP of the methane–water dimer (Figure 4.15b) because the polarity of the water H atom dominates the MEP in that region. But if we remove the water molecule, the net increase in negative charge in the methane molecule where it interacts with the partially positive H atom of the water is apparent (Figure 4.16).

Making the figure
This is a not a plot one makes every day, so the process is a bit involved. The main trick is to construct the methane part of the density from the corresponding localized molecular orbitals, and then tricking GAMESS into printing a file with just those LMOs that MacMolPlt can read. The procedure is described in some detail Section 5.5 of the book, so here I just post the corresponding screencast.

Canonical and localized molecular orbitals

Figure 4.10. 0.045 au isosurfaces of the four valence canonical MOs of NH₃ computed using B3LYP/6-31G(d). A value of 0.045 au is chosen because it results in a 0.002 au isodensity surface when squared.
From Molecular Modeling Basics CRC Press, May 2010.
Click on the picture for an interactive version

According to Valence Shell Electron Pair Repulsion (VSEPR) theory many bond angles involving elements such as C, N, and O are close to 109.5^o because the four valence electron pairs that surround these atoms adopt a tetrahedral geometry to minimize repulsion. In the case of CH₄ the H-C-H angle is exactly 109.5^o because the repulsion between the four electron pairs in the C-H bonds are identical. NH₃ has a lone pair that is fatter than a bond 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 of 107^o, slightly smaller than 109.5^o.

However, this is far from obvious when looking at the four valence MOs of NH₃ (Figure 4.10) computed using B3LYP/6-31G(d). The reason is that that MOs are not unique, and that these MOs (which lead to a diagonal Fock matrix and are known as canonical MOs) are not the MOs where the electron pair repulsion is a minimum. Algorithms have been implemented that find a new set of MOs (localized MOs or LMOs) that represent a linear combination of canonical MOs for which the MO–MO repulsion is a minimum (Figure 4.11).

Figure 4.11. 0.045 au isosurfaces of the four valence localized MOs of NH₃ computed using B3LYP/6-31G(d). There are three N–H bond LMOs [(a)–(c)] and one lone pair LMO [(d)].
From Molecular Modeling Basics CRC Press, May 2010.
Click on the picture for an interactive version

LMOs for which the inter-orbital repulsion is a minimum are called energy localized orbitals or Edmiston-Ruedenberg orbitals. Other popular choices are Foster-Boys and Pipek-Mezey LMOs, which use different localization criteria.

Using MacMolPlt
The screencast below shows how to compute Ruedenberg LMOs for ammonia using GAMESS (local=ruednbrg in the $contrl group), and how to display the LMOs, as well as the canonical MOs, in MacMolPlt. I also show how to identify the HOMO and LUMO canonical MOs in MacMolPlt.

Note that I specify a geometry optimization in the GAMESS input file. GAMESS will only compute the LMOs for the optimized geometry.

Using Jmol
I use Jmol for the interactive figures and the scripts can be found here and here. Unlike the density and electrostatic potential, Jmol can generate its own grid data, so to display MO number 2 you simply use "mo 2" in the script. I use

mo 2; mo cutoff 0.045; mo fill nomesh; mo translucent 0.2

to make it a little prettier.

Jmol cannot find the MOs in a GAMESS geometry optimization file, and Jmol only stores the LMOs if present. So to display the canonical MOs with Jmol you need a single point energy calculation (output file) and to display the LMOs with Jmol you need a single point energy calculation with local=ruednberg added (output file).

Common error messages in GAMESS: Failure to locate stationary point

***** FAILURE TO LOCATE STATIONARY POINT, TOO MANY STEPS TAKEN *****
  UPDATED HESSIAN, GEOMETRY, AND VECTORS WILL BE PUNCHED FOR RESTART
**** THE GEOMETRY SEARCH IS NOT CONVERGED! ****

This error message was produced by the following input file. Can you see what's wrong?

$contrl runtyp=optimize icharg=1 $end
 $basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

Actually, there is no problem with the input file as such. A geometry optimization is an iterative process and if the gradient it not below the convergence criteria within 20 steps, GAMESS will stop and print out the message shown above.

The solution is simply to take the last set of coordinates and run the optimization again, as I show in the screencast below.

As I've mentioned in a previous post I think the default criterion for geometry convergence (0.0001) is too strict, and the default number of steps (20) is too small. So I usually use 0.0005 and 50, respectively.

$statpt nstep=50 opttol=0.0005 $end 

Common error messages in GAMESS: charge and multiplicity

  *** CHECK YOUR INPUT CHARGE AND MULTIPLICITY ***
THERE ARE    33 ELECTRONS, WITH CHARGE ICHARG=  0
BUT YOU SELECTED MULTIPLICITY MULT=  1

This error message was produced by the following input file. Can you see what's wrong?

 $contrl runtyp=optimize $end
 $basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

As GAMESS tells you, it is not possible to have singlet state (i.e. a multiplicity of 1) if you have an odd number of electrons. You need to tell GAMESS whether you want to change the multiplicity (for example, $contrl mult=2 scftyp=rohf $end) or add ($contrl icharg=-1 $end) or remove ($contrl icharg=1 $end).

When looking at the molecule it is pretty clear that you want to remove an electron, since ammonium has a positive charge. So the first line in the input should be changed to

 $contrl runtyp=optimize icharg=1 $end

Click on the picture for an interactive version

Notice that you don't have to specify that the positive charge is on ammonium. The quantum mechanics will take care of that.

Common error messages in GAMESS: error reading IDUM

 **** ERROR READING VARIABLE IDUM     CHECK COLUMN  1
H     1.0    -0.08386     2.76975     0.70945                                
....V....1....V....2....V....3....V....4....V....5....V....6....V....7....V....8

This error message was produced by the following input file. Can you see what's wrong?

 $contrl runtyp=optimize $end
$basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

The answer is the missing space in front of the $basis, and the solution is to add a space in front of it, just like for $contrl and $basis.

Why does GAMESS complain about a line in $data when the problem is with $basis? The reason is that it is possible to specify a basis set for each atom in $data, so GAMESS tries to read this information if no $basis group is found.

The error message would clearly benefit from the line: "Did you mean to specify a $basis group?"

Common error messages in GAMESS: No $data

**** ERROR, NO  $DATA   GROUP WAS FOUND
EXECUTION OF GAMESS TERMINATED -ABNORMALLY-

This error message was produced by the following input file. Can you see what's wrong?

 $contrl runtyp=optimize $end
 $basis gbasis=pm3 $end
$data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

The answer is the missing space in front of the $data, and the solution is to add a space in front of it, just like for $contrl and $basis.

Electron density

Figure 3.1. (a) A contour plot of the density of H₂O computed using RHF/6-31G(d). The maximum and minimum contour values are 0.05 and 0.002 aus. (b) The corresponding 0.002 au isodensity surface. (c) The surface corresponding defined by atomic spheres with van der Waals radii.
From Molecular Modeling Basics CRC Press, May 2010.

Here's a screencast about how I made the figure:

When making the contour plot (Figure 3.1.a) I pick 0.05 au as the maximum value (this will be the contour line closest to the nuclei) and 25 contour lines. This means that the spacing between the contour lines and thus the outer contour line will be 0.05/25 = 0.002 au.

Another, more common, representation of the density is a 3D version of one of the contour values: the isodensity surface (Figure 3.1.b). A common choice is 0.002 au, since that corresponds roughly to experimental estimates of molecular size, such as the van der Waals surface (Figure 3.1.c).

Unfortunately, MacMolPlt doesn't have a van der Waals display style, so I have to use Avogadro. This means I have to re-size (by eyeball) this part of the figure to make it the same size as the density plots. See this post about using screen capture to get a file with the graphic.

The interactive version of this figure is the subject of a future post.

The trouble with transition metals II: SCF convergence

In a previous post I showed how to build two structural models containing transition metals: a small model of a zinc finger and a ferrocene. This post is about computing the energy using GAMESS.

The first challenge in computing the energy is choosing the correct charge and multiplicity. Most organic molecules are singlets (i.e. all orbitals are doubly occupied) and have an easily identifiable charge. This is often not the case with molecules containing transition metals. Many transition metal atoms have unpaired d electrons and can have several different charges (oxidation states). Also, the charge of the ligands is not always obvious, though they almost always are singlets.

The zinc finger model is relatively easy: zinc is almost always Zn²⁺ and has 5 doubly occupied d orbitals, so that's a pretty safe bet. The zinc finger model contains two neutral groups (the imidazole rings) and two negatively charged groups (CH3-S^-). So the charge of the entire system is zero and the multiplicity is 1.

The ferrocene model is a bit more complicated: iron is commonly found in both the +2 and +3 oxidation state. Before one starts the calculation it is important to find out which one is appropriate for ferrocene. Google says +2 (meaning 6 d electrons) and all doubly occupied orbitals (see the orbitals at the bottom of the page), so a singlet. Furthermore, it is also important to know that each ring has a -1 charge, so the total charge of the complex is 0 and the mutiplicity is 1.

Now, I was all set to talk about SCF convergence problems, i.e. that the zinc finger would converge fine but the ferrocene would give problems. But when I ran a RHF/6-31G(d) energy calculation

$contrl runtyp=energy $end
$basis gbasis=n31 ngauss=6 ndfunc=1 $end
$scf dirscf=.t. $end

both runs converged fine!

"Luckily" when using a smaller basis set (RHF/STO-3G) leads to convergence problems - for both. Here is the relevant part of the SCF output for zinc finger.

                                                                                                                                                                     NONZERO     BLOCKS
ITER EX DEM     TOTAL ENERGY        E CHANGE  DENSITY CHANGE     ORB. GRAD      INTEGRALS    SKIPPED
 1  0  0    -3064.6703808712 -3064.6703808712   1.733010527   0.000000000        9751402     457807
 2  1  0    -3059.0602240553     5.6101568159  16.724847696   0.794728909        9737106     471195
 3  2  0    -3053.0258026844     6.0344213709  16.846293078   0.656307634        9778274     459644
 4  3  0    -2935.9005712044   117.1252314800 115.104028436   1.822663152        9792060     446482
 5  4  0    -2923.8944526114    12.0061185930 115.740577390   0.620923864        9792989     441267
 6  5  0    -2768.2361993103   155.6582533011 115.747362137   5.442595435        9783876     443565
 7  6  0    -2921.1390148151  -152.9028155048 115.772970667   0.810505145        9782321     444765
 8  7  0    -2775.0326048138   146.1064100013 115.770983819   4.741230605        9783214     444582
 9  0  0    -2920.9171138162  -145.8845090024  60.570511833   0.796143103        9783475     444487
10  1  0    -3003.9501726068   -83.0330587906  54.636432025   0.749401761        9790866     444493
11  2  0    -2862.1853212282   141.7648513786 115.095938299   1.824897836        9796139     439522
12  3  0    -2921.4368542174   -59.2515329892 115.495602218   0.722894426        9783664     443903
13  4  0    -2848.2964423996    73.1404118178 115.538938253   2.394847457        9782038     444795
14  5  0    -2905.4155965583   -57.1191541586 115.529755578   0.678236838        9781255     445351
15  0  0    -2859.1068450197    46.3087515386  14.660262004   1.964358011        9783437     444643
16  1  0    -3019.7713654862  -160.6645204665  14.032213984   0.801579670        9794434     443563
17  2  0    -2934.9445251088    84.8268403773  99.891887745   1.748565491        9794346     445610
18  3  0    -2905.6599083961    29.2846167127 100.392269592   1.167341691        9789953     445157
19  4  0    -2864.9366576173    40.7232507788 104.001625891   1.797814244        9781805     446591
20  5  0    -2855.1081148594     9.8285427580 103.862237046   1.447532031        9781260     446026
21  6  0    -2865.3401581051   -10.2320432457  82.200569407   1.143844678        9780463     446052
22  0  0    -2857.7261418857     7.6140162194  78.128990865   1.766124007        9781544     445447
23  1  0    -3036.7926259046  -179.0664840190   5.185277571   0.905312779        9796083     442344
24  2  0    -2977.1610659805    59.6315599241   5.419478057   1.193389134        9796601     444172
25  3  0    -2859.5939849458   117.5670810347 114.799971576   1.724872824        9789000     444821
26  4  0    -2913.7281172708   -54.1341323250 115.553208606   1.094842493        9785096     443663
27  0  0    -2845.9860251842    67.7420920865  16.890315601   1.995942532        9783715     444171
28  1  0    -2989.2091196251  -143.2230944409  16.260193167   1.244487455        9796844     439846
29  2  0    -2911.6240721064    77.5850475187 114.870828357   1.711669748        9798193     439590
30  3  0    -2897.8498052725    13.7742668339 115.443658127   0.704588469        9789911     442464

SCF IS UNCONVERGED, TOO MANY ITERATIONS
  TIME TO FORM FOCK OPERATORS=      84.7 SECONDS (       2.8 SEC/ITER)
  FOCK TIME ON FIRST ITERATION=       2.6, LAST ITERATION=       3.0
  TIME TO SOLVE SCF EQUATIONS=       0.4 SECONDS (       0.0 SEC/ITER)

FINAL RHF ENERGY IS        0.0000000000 AFTER  30 ITERATIONS

Clearly, the this SCF is not converging and giving it more iterations is not going to make the problem go away. Instead, different algorithms for SCF convergence are needed and

$scf dirscf=.t. diis=.t. damp=.t. $end

leads to convergence in for both ferrocene and zinc finger.

DIIS stands to Direct Inversion of the Iterative Subspace, which is an extrapolation procedure, while DAMP "damps" any oscillations between in the SCF.

There are a few other tricks to SCF convergence but they are clearly not needed here. If DIIS and DAMP doesn't solve the problem for you, leave a message describing your molecule and I'll see what I can do.

Get a reaction: Intrinsic Reaction Coordinate

Lucas requested a post on computing an intrinsic reaction coordinate (IRC). An IRC is a special case of the minimum energy path (MEP) connecting reactant, transition state, and product (more about this at the end of the post).

There are many uses for the IRC, the most basic of which is to verify that the transition state you found actually connects the reactants and products you think it does (this is not always obvious from the normal mode associated with the imaginary frequency). It is also a good way to identify any reaction intermediates you haven't thought about. The IRC is also a way generate an animation of the reaction (see for example ChemTube3D).

An IRC is generated by following the gradient down-hill from the transition state to products and reactants. (It's a bit like a geometry optimization initiated at the transition state, but the algorithm used to follow the gradient is different). Computing a complete IRC thus requires a minimum of two separate runs, and you must find the TS before you can compute an IRC.

The screencast shows how to generate the IRC associated with the 3TSa transition state I found in a previous post. The first input file I generated contain the 3TSa coordinates and the following keywords:

$contrl runtyp=irc $end
$basis gbasis=pm3 $end
$irc saddle=.t. tsengy=.t. forwrd=.t. npoint=100 stride=0.2 opttol=0.0005 $end

saddle=.t. tells GAMESS that the coordinates in $data corresponds to a TS. Because the TS is a stationary point the gradient is zero (or very small), so there is no gradient to follow down-hill. Therefore, the IRC is initiated by displacing the geometry along the normal mode corresponding to the imaginary frequency. Therefore, saddle=.t. GAMESS expects a $HESS group, which contains the frequency information.

Whether the normal mode points towards the reactants or products is completely arbitrary, and the forwrd=.f. keyword is used to change the direction in which the normal mode points. forwrd=.t. is the default and I only include it in the first run for completeness.

tsengy=.t. tells GAMESS to re-compute the energy of the TS (re-compute because you already computed it when you computed the frequencies). If you want to make an energy plot in MacMolPlt as I show in the screen cast, you must set tsengy=.t.

The displaced geometry will have a non-zero gradient, the negative of which points downhill. The geometry is displaced in this direction, the gradient is recomputed, and the procedure is repeated npoint=100 times or until the (root-mean-square) of the gradient is below opttol=0.0005 atomic units, at which point you are at a minimum. (Notice that the default for npoint is 1, and must always be changed).

stride=0.2 tells GAMESS how far along the gradient it should displace at each step (in this case 0.2 atomic units). The default is 0.3, but when I tried this, the IRC calculation in the forward direction failed after 7 steps, because the coordinate deviated too much from the minimum energy path, so I had to decrease it. This means I will have to take more steps to reach the minimum.

As you saw in the screencast, GAMESS takes 100 steps in each direction, which means the IRC does not go all the way to the minimum on each side of the TS. GAMESS does print out restart information to continue the IRC, which I didn't do in this case.

Finally a note on MacMolPlt. As I show in the screencast, MacMolPlt can be used to combine the two IRC runs (forwrd=.t. and forwrd=.f.) by using Files > Add Frames from File... You should first load the file that goes to product and then add the file from the IRC that goes to reactants. Here you should tell MacMolPlt to reverse the order of the structures in the file, so the animation goes from reactions to transition state to product.

This is done by selecting "make these points negative". "Negative" refers to the x-coordinate in the energy plot. The x-coordinate is what separates an IRC from a minimum energy path. An IRC is a plot of the energy of each structure along the minimum energy path vs the root-mean-square change in the mass weighted Cartesian coordinates relative to the TS structure. The TS thus has a x-coordinate of zero, while structures leading to the reactants are defined by negative RMSD values.

If the energy is plotted against something else, such as a bond-length, then it is simply a minimum energy path. Because an IRC uses mass-weighted coordinates, using different isotopes will lead to different IRCs, and one use of IRC is to determine tunneling corrections to the reaction rate.

Get a half-life: From transition state to rate constant

[source]

In a previous posts I showed how to find the TS for hydrolysis of 3 at the RHF/3-21G level of theory. In this post I show how to compare the computed results with the experimentally observed half-life (t_1/2) of 88 hours.

The half-life is connected to the the rate constant for the reaction (k) by

(1)

The rate constant can be estimated from the activation free energy (ΔG^≠) by transition state theory:

(2)

(here the units are those of a unimolecular reaction). The activation free energy is the free energy of the TS relative to the reactant

ΔG^≠ = G(3TS) - G(3) = ΔE^≠+ ΔG^≠_TRV

Which can be rewritten as the difference in electronic energy and the free energy correction due to translation, rotation, and vibration. I showed how to extract the RHF/3-21G electronic energy (-413.5161 Hartrees) and free energy correction (81.841 kcal/mol) for the 3TSa in a previous post.

In the screencast at the end of the post I show how to obtain the corresponding values (-413.6014 Hartrees and 82.889 kcal/mol) for 3. So,

ΔG^≠(3TSa) = (-413.5161 - (-413.6014))*627.51 + (81.841 - 82.889)
= 52.5 kcal/mol

You can get the corresponding experimental number from the half life (88 hours) and equations (1) and (2): 25.1 kcal/mol and, yep, the RHF/3-21G activation free energy is way off. The corresponding barrier for 3Tsb (51.0 kcal/mol) is only a little better.

There could be many reasons why the activation free energy is so severely overestimated. Two of the most likely suspects are:

1) Level of theory. TS structures, and therefore the underlying wavefunctions, tend to be more complicated than the minima structures. Thus, larger basis sets and electron correlation methods may be needed to obtain accurate activation free energies. This tends to lower the energy of the TS more than the minima, thus lowering the barrier compared to lower levels, but this is not always the case.

2) A different reaction mechanism. What might seem the most straightforward path from reactants to products is not always the one with the lowest barrier. For example, in the case of amide hydrolysis a stepwise mechanism is also possible in which the hydroxyl hydrogen is transferred to the carbonyl O atom to create a (chiral) intermediate.

TSs corresponding to such a mechanism must also be found to see if they lead to lower activation free energies. If you followed the last few posts of finding TSs you now have all the skills needed to do that.


Here is the screencast that shows how to find the energies for 3:

I paste in the following set of keywords

$contrl runtyp=optimize $end
$system mwords=125 $end
$statpt opttol=0.0005 nstep=50 hssend=.t. $end
$force nvib=2 $end
$contrl nzvar=1 $end
$zmat dlc=.t. auto=.t. $end

In the RHF/3-21G optimization + frequency run I forgot to add

$scf dirscf=.t. $end

which would have speeded up the calculation, as I described in a previous post.

GAMESS, memory and parallel

The screencast above is a repeat from the last post, where I computed the frequencies of a molecule at the RHF/3-21G level of theory (second screencast) and discovered that the amount of memory that GAMESS requests (1,000,000 words) was not enough. While GAMESS tells me how much it needs in the output file, it does so after a lot of computation, which is wasted because I have to start over.

In this post I discuss the basics of GAMESS memory requirements. The simplest case is when you are not running in parallel, i.e. the number of processors in GAMESSQ is set to 1, so I discuss this case first.

1. Memory is specified in $system, and the default is 1 mega-word (i.e. 1,000,000 words):

$system mwords=1 $end

mword can only be integers (1, 2, 3, ..). (In the screencast above I used the older "memory" keyword but mwords is much more convenient.)

2. A word is 8 bytes, so the default memory request (1 mega-word) is a very modest 8 MB of RAM.

3. This is the maximum amount of memory GAMESS is allowed to use. Depending on the type of calculation GAMESS might use less. The option is there to avoid filling up the memory on the computer entirely, which will crash or freeze the computer.

4. Here is the simplest way to deal with memory: My current laptop computer has 2 GB of RAM, and I use it for other things while GAMESS is running, so I am willing to give GAMESS a maximum of roughly 1 GB of RAM. This translates to

1 GB ≈ 1,000 MB ≈ 125 mwords

5. Adding

$system mwords=125 $end
to all input files will allow me to run most GAMESS jobs that I would want to run on a laptop without ever getting a "not enough memory" error (of course if you have less memory you need to adjust mwords accordingly).

That's basically it. What follows is some more details that most casual users of GAMESS won't have to worry about.

6. Most common GAMESS runs will never get close to using 1 GB of RAM, meaning the memory is free to be used by other programs. The most common types of runs that potentially could use that much memory are frequency calculations using RHF and any kind of MP2 calculation. Slightly less common ones are TDDFT and NMR calculations. In GAMESS, DFT frequency calculations are done numerically, and do not require a lot of memory.

7. If you are in doubt whether you have reqeusted enough memory to perform a calculation it is possible to use GAMESS to check using

$contrl exetyp=check $end

This keyword will simulate an actual GAMESS run by skipping the most time consuming steps so it is very fast.

I show two examples (an RHF/3-21G frequency and MP2/6-31G(d) single calculations on 3TSa) in this screencast.

Parallel runs

8. Many desktop computers (and almost all larger computers) have more than one CPU (also known as cores). For example, my current laptop has one processor with 2 cores. Thus, I could make my GAMESS calculations go faster by specifying 2 processors when submitting with GAMESSQ (processor = core in GAMESS-speak). This means that two separate but related GAMESS calculations are running simultaneously, and this affects the memory request:

9. Most types of runs use replicated memory. This means that if mwords=1 but I ask for 2 processors, GAMESS will use a maximum of 2 mega-words. Thus, if I routinely use 2 processors when running on my laptop and want to impose a limit of 1 GB RAM, I should use

$system mwords=63 $end

instead of 125.

10. The most common exception to this is MP2 calculations, where GAMESS uses both replicated (mwords) and distributed memory (memddi). Distributed memory is memory shared by the cores. If I specify memddi=100 and ask for 2 processors then 100 mega-words of RAM is distributed among the 2 processors (the simplest case is that each core gets 50 mega-words, but GAMESS figures that out for you).

So running on 2 processors with

$system mwords=15 memddi=100 $end
requests a total of (15 + 15 + 100 =) 130 mwords. You can use exetyp=check to figure out the optimum values of memory and memddi, as I show in this screencast.

The screencast shows how you can use 1 processor and

$contrl exetyp=check $end
$system memddi=100 parall=.t. $end

to check the memddi requirements. You need to make memddi large for this to work, but GAMESS will not use this memory during the check-run. Note that parall only "does something" for check runs: true parallel runs are specified by choosing more than 1 processors when you submit GAMESS.

You can get a complete list of runs that require memddi in Chapter 2 of the GAMESS documentation under the entry for $system. Also, Chapter 5 has a more in-depth description of how memory is handled.

Finding a transition state: amide hydrolysis 2

In a previous post I found two TSs, called 3TSa and 3TSb, at the PM3 level of theory.

In this post I show how the PM3 structure and frequency information can be used as a starting point for finding 3TSa at the RHF/3-21G level of theory.

I use RHF/3-21G as an example of an expensive ab initio method where we want to keep the number of frequency calculation to a minimum, since it takes a long time to compute.

So I extract the PM3 optimized geometry from the output file of the TS search, and add to this the corresponding frequency information in the form of the $HESS group from the .dat file.

The keywords I paste in look like this (many of these keywords have been described in a previous post):

$contrl runtyp=sadpoint $end
$scf dirscf=.t. $end
$statpt opttol=0.0005 nstep=50 $end
$statpt hess=read $end
$contrl nzvar=1 $end
$zmat dlc=.t. auto=.t. $end

Note the $basis group is missing, but that I specify it using Avogadro. It's important to pick the basis set before pasting the remaining keywords in because Avogadro deletes some of the lines in the file when you select something in the menu.

Based on the PM3 guess, GAMESS finds the RHF/3-21G TS in 29 steps.

Actually, I can't be 100% sure the structure is a TS because I haven't verified that the structure has an imaginary frequency. Unlike with PM3 I don't automatically compute the frequencies at the end of the TS search (using hssend=.t.).

The reason is that it is a good idea to verify that the structure looks like a TS before committing CPU time to an expensive frequency calculations. As I showed in a previous post it can happen that TS searches end up find the reactants and products.

Luckily the TS structure has one, and only one imaginary frequency as you can see in this screencast

As you saw the default amount of memory that GAMESS requests (1,000,000 words) is insufficient for the RHF/3-21G frequency calculation, but it tells you how much it needs (~1,100,000 words). A "word" is 8 bytes, so we are talking about 10 MB here. I will discuss this issue in more details in a future post.

Finally, in principle it is possible to adjust the memory keyword in Avogadro (in the Advanced setup), but that option is currently not working.

Now we need to repeat this procedure for 3TSb. The procedure are exactly the same, so I haven't made a screencast. The TS search takes 36 steps.

GAMESSQ for everyone

When I made the post about GAMESSQ, I had no idea it was only being distributed for the Mac (together with the Mac binary) - I usually don't make posts on software that is not cross-platform.

However, now Brett Bode has made a general release of GAMESSQ. Here you can find a Windows executable and the source code, the later being the only current option for Linux users.

It is not inconceivable that one or two Linux binaries will be available soon. What flavor if Linux would be most useful to MolModBasic readers?

Finding a transition state: amide hydrolysis I

[source]

In a previous post I showed how to find a TS for an Sn2 reaction involving a small molecule (it's a good idea to look at this post now if you haven't seen it before, and check out this post as well).

In this and future posts I discuss the steps involved in finding a TS for a reaction involving complicated molecules shown in the figure above (I have shown how to build some of the reactants in a previous post).

Say you want to find the TS for the hydrolysis of 4 and 5 at a relatively high level of theory. My advice is to start by finding the TS at a low level of theory (PM3) for a smaller representative reaction (hydrolysis of 3). This is the subject of today's post.

The basic procedure is the same as in the previous post:

1. I build a TS-like structure by constraining a key distance during an energy minimization.
2. I check whether the optimized geometry has an imaginary frequency that looks like the reaction coordinate
3. If it does, I use it as a starting point for the TS search.

In the previous example, one bond broke as the other one formed, so I simply forced one bond (the C-F bond) to form by constraining it to 2.0 Å and the optimization took care of the rest.

In the hydrolysis of 3, two bonds have to break (the C-N and O-H bonds) as two other bonds form (the C-O and N-H bonds), so what to constrain? In general I have found it a good idea to try to construct an intermediate structure, such as "3TSa" shown below.

Here the O-H bond is already broken and the C-O and N-H bonds are already formed. The question is now reduced to whether I should try to break the C-O or the C-N bond. I opted for the C-O bond which I start breaking by constraining it to 2.0 Å. The reason for choosing the C-O bond is that I assumed that the (potential) energy of 3 is lower than the products (because of ring strain) so the TS should be more product-like according to Hammond's postulate (you can of course check this computationally).

Why do I pick 2.0 Å? Well, you have to start somewhere and this is a nice round number (for bonds involving 3rd row atoms I would use 2.5 Å). As you see in the screencast below, 2.0 Å doesn't work but the TS search tells me how to change it, i.e. increase it to make it more reactant like as the TS search goes to products.

As you saw in the screencast I have pre-typed some text with keywords that I use again and again, and simply paste the text in the input file and alter it when necessary. For the constrained optimization the text looks like this:

  $contrl runtyp=optimize $end
  $basis gbasis=pm3 $end
  $statpt opttol=0.0005 nstep=50 hssend=.t. $end
  $force nvib=2 $end
  $contrl nzvar=1 $end
  $zmat dlc=.t. auto=.t. ifzmat(1)=1,x,x fvalue(1)=2.00 $end

and for the TS search the text looks like this.

  $contrl runtyp=sadpoint $end
  $basis gbasis=pm3 $end
  $statpt opttol=0.0005 nstep=50 $end
  $statpt hess=calc ihrep=5 hssend=.t. $end
  $force nvib=2 $end
  $contrl nzvar=1 $end
  $zmat dlc=.t. auto=.t. $end

You might also notice that I am being a "CPU slob" in the sense that I redo some computations I've already done. For example, I compute the frequencies for the optimized constrained geometry twice: at the end of the optimization and at the beginning of the TS search. That is because I am using PM3, which is such an efficient method that the computer can re-compute the freqeuncies in less time than it takes me to transfer the $HESS group to the TS search input file (as shown in a previous screencast)! Similarly, when I find that 50 steps is insufficient and just restart with nstep=100 steps, which recomputes the first 50 steps. Good thing sloth is not a deadly sin.

Another advantage of constructing a hypothetical intermediate in this case is that you realize (hopefully) that the TS contains a chiral center, even though the reactants and products are achiral. So you have to look for two TSs, to see which one is lowest in energy.

The next screencast shows how I find 3TSb, which is made easier by the fact that I have a good idea of what the constraint value should be (2.32 Å) from the 3TSa geometry I found previously.

Finding a transition state: an Sn2 reaction

Here is a screencast I made about finding a transition state for the reaction

F^- + CH₃Cl -> Cl^- + CH₃F

The general approach is as follows

1. Build the structure (F^- + CH₃Cl)
2. Optimize the geometry while constraining a key distance (between the F and C atom) to a certain value (2.0 Å)
3. Compute the frequencies for the optimized geometry to check for an imaginary frequency
4. If there is an imaginary frequency, use this geometry and frequency information as a starting point for the TS search

The distance constraint is imposed using

$zmat ifzmat(1)=1,1,6 fvalue(1)=2.0 $end

which is used in conjunction with the delocalized internal coordinates I discussed previously. 1,1,6 defines the distance (denoted with the first "1") between atoms 1 and 6. Note that this distance does not have be 2.0 Å in the guess geometry (though the closer it is to 2.0 Å the faster the optimization will converge).

I also show how the frequency information contained in the $HESS group, found at the end of the .dat file that GAMESS produces along with the output file, can be transferred to the input file for the TS search, where hess=read is used in $statpt. Any text editor can be used for this. Note that the hssend=.t. keyword has to be used to cause the frequencies to be computed at the end of the optimization. I have showed the use of this keyword in a screencast in a previous post.

The complete $statpt groups looks like this

$statpt opttol=0.0005 nstep=50 hess=read ihrep=5 hssend=.t. $end

ihrep=5 causes the frequencies to be recomputed at every 5th step. This really helps convergence but is computationally expensive.

Note that there are no constraints imposed when searching for the TS.

The key question is of course why I chose the C-F distance and 2.0 Å for the constraint value. This is simply the first thing that popped into my head and it happened to work. If it hadn't, I would have tried some other values and, if none worked, I would have tried to the C-Cl distance. Future posts will show examples of trying difference distances.

Because you may have to try different constraints before succeeding, it is a good idea to find the TS at a cheap level of theory, such as PM3. The C-F distance in the PM3 TS structure (2.10 Å) should provide a good guess for TS searches at higher levels of theory.

I tried it for this reaction using RHF/3-21G, but with no luck. That is because I forgot the very first step in searching for a TS: make sure there is one! In the screencast below I show how a RHF/3-21G geometry optimization (with no constraints) starting from F^- + CH₃Cl leads to Cl^- + CH₃F, meaning that the former structure is not a minimum on the gas phase RHF/3-21G potential energy surface. So there is not much point in looking for the corresponding TS.