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.

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).

I'm positive

Figure 4.2. 0.002 au isodensity surface with superimposed electrostatic potential of (a) Li⁺, (b) Na⁺, and (c) K⁺ ion. The maximum potential value is 0.8 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.

Here is an example of how computational chemistry can be used to enhance teaching at the general chemistry level.

Why does the ionization energy decrease on from Li to Na to K? That's the same as asking why the electron affinities decrease on going from Li⁺ to Na⁺ to K⁺.

Figure 4.2 show the ions colored by how positive they are at the surface (i.e. the electrostatic potential superimposed on the 0.002 isodensity surface). The darker the
color the more positive the ion, and it is clear that the Li+ ion is “more positive” than Na⁺, which is more positive than K⁺. Thus, more energy should be released when adding an electron to Li⁺ compared to Na⁺, and hence more energy is needed to remove an electron from Li
compared to Na (and similarly for K).

The reason why Li⁺ is “more positive” than Na⁺ or, more accurately, why the potential on the 0.002 au isodensity surface is more positive for Li⁺ than for Na⁺, is that the former is a smaller ion than the latter, so the surface is closer to the +1 charge at the center of the ion.

Figure 4.4. 0.002 au isodensity surface with superimposed electrostatic potential of (a) Li⁺, (b) Ne⁺, and (c) Na⁺ ion. The maximum potential value is 0.8 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.

This rationalization is, of course, only qualitative and is is not predictive of the ionization energies among different groups. For example, based on Figure 4.4 one would expect that Ne would have roughly the same ionization potential as Na, which is not true at all.

When making these figures it is very important to get the relative sizes of the ions correct, but this can be difficult since each image can be zoomed to an arbitrary size. The screencast below shows how to control this in MacMolPlt using the Manual Windows Parameter window.

The same point can also be made with a "quick-and-dirty" electrostatic potential, as shown in this interactive figure. Here the electrostatic potential is due to a plus one charge centered at the atom (rather than the nuclear charge and the electron density) and the surfaces are the spheres defined by empirical ionic radii. The Jmol script can be found here, and the mol2 files here, here, and here.

Click on the picture for an interactive version

I have my moments

Figure 3.7. Contour plot of the RHF/6-31G(d) electrostatic potential and 0.002 au isodensity surface of (a) CH₃COO^-, (b) HF, and (c) F₂. The maximum/minimum contour values are, respectively, 0.5/0.025; 0.1/0.005; and 0.005/0.00025 au respectively. Blue corresponds to a negative potential. In each case the outer-most contour looks like the corresponding contour in the electrostatic potential due to a charge, dipole, and quadrupole, respectively.
From Molecular Modeling Basics CRC Press, May 2010.

This figure makes 3 points:

1. It shows what the electrostatic potential of a charge, dipole, and quadrupole looks like.

2. It show the relative strengths of the electrostatic potentials due to a charge, dipole, and quadrupole.

3. It shows that the electrostatic potential of a charge, dipole, and quadrupole deviates significantly from the actual electrostatic potential near the molecular surface.

Here is a screencast showing how I made Figure 3.7a. (If I were to do it over I would have chosen red for negative and blue for positive... oh, well.)

Here is an interactive version of the figure. Click on the picture to load it. Remember it's Jmol so you can rotate it and zoom as you like (Mac users: this works best with Safari).

Figure 3.7. The RHF/6-31G(d) electrostatic potentials of acetate, HF and F₂.
Click on the picture for an interactive version

See these two posts (here and here) on how to make plots like that with Jmol. From a Jmol perspective the only new thing is that I show two isosurfaces simulateneously. This is done by loading the file twice to create two "frames" that Jmol can display simultaneously. You can find the script file with all the commands here, but the general syntax is:

load files file1.xyz file2.xyz
frame 1.1; isosurface surf1 plane {0 0 0 0} contour 40 color range -0.05 0.05 "potential.cube.gz"
frame 2.1; isosurface surf2 0.002 "density.cube.gz" 

Even though you want a 2D contour plot of the potential, it is necessary to make a 3D cube file. Because I use unusually low cutoffs to show the outer contours, I had to trick MacMolPlt into making a bigger grid. I show how on the screencast below.

Electrostatic potential maps

Figure 3.4. (a) Contour plot of the electrostatic potential of H₂O. The maximum and minimum contour value is 0.5 Hartrees/electronic charge (au). (b) The corresponding 0.05 au isopotential surface. (c) The electrostatic potential displayed on the 0.002 au isodensity surface of water. The maximum (darkest blue) value corresponds to 0.05 au.
From Molecular Modeling Basics CRC Press, May 2010.

Here is a screencast of how I made the figure:

The files I used were created in a previous post. The various values specified above are determined using trial and error, i.e. I kept fiddling with it until it looked good to me. When using plots like this it is important to specify these values, because using different values can lead to very different looking plots.

Also, different programs use different color scales and color intensity to indicate positive and negative charge, so it is rarely possible to directly compare electrostatic potential maps from different programs.

See this post about using screen capture to get a file with the graphic. This was done for each plot, and the combined in a word processor.

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

Electron density reloaded

Here is an interactive version of a figure I described in a previous post. Click on the picture to load it. Remember it's Jmol so you can rotate it and zoom as you like (Mac users: this works best with Safari).

Figure 3.1. The RHF/6-31G(d) electron density of water.
Click on the picture for an interactive version

The Jmol animation loads this file (h2oprinc.xyz), which looks like this

3
jmolscript: script "http://propka.ki.ku.dk/~jhjensen/h2odensity.spt"
O -0.0000 0.0643 0.0000
H -0.7541 -0.5091 0.0000
H 0.7541 -0.5091 0.0000

and which, in turn, loads a script (h2odensity.spt), which looks like this

isosurface planex plane {0 0 0 0} contour 20 color absolute 0.002 0.05 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
delay 3
spin y 20
delay 10
spin off
isosurface planey plane {1 0 0 0} contour 20 color absolute 0.002 0.05 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
spin y 20
delay 10
spin off
isosurface planez plane {0 1 0 0} contour 20 color absolute 0.002 0.05 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
spin x 40
delay 10
spin off
isosurface threed 0.002 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
color isosurface red ; color isosurface translucent 0.15
spin y 20
delay 10
spin off
color isosurface red ; color isosurface translucent 0.5
select all; spacefill 100 %babel
spin y 20
delay 10
spin off

The screencast below shows how I made the cube file that contains the electron density information using MacMolPlt. The program I used to convert the MacMolPlt file to a cube was written by Jonathan Gutow and can be download here. In the screencast I mistakenly named the cube file h2odensity.cube.gz, but that's easy to change.

I start by loading the RHF/6-31G(d) optimized geometry I computed for a previous post, and reorienting. The orientation makes it easier to define the plotting plane for the contour plots.

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.

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.

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.

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.

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.

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.

Tools of the trade

Here are some wonderful programs that I use all the time in teaching and research (all programs are free and work on Windows, Mac, and LINUX):

I use Avogadro for building molecular structures and setting up input files for GAMESS (see below). Avogadro is still in beta version but quite close to "1.0" release. The main practical result of this is that there is very little documentation, and no real manual. Fortunately, Avogadro is very intuitive and there are a few screencasts that are quite informative (though I say so myself - some of them are mine).

I use GAMESS for quantum chemical calculations. It is not yet possible to submit GAMESS input files created in Avogadro automatically to the GAMESS program, so you must save the input file in Avogadro and submit the GAMESS job from the command line.

I visualize the GAMESS results by opening the output files in MacMolPlt (which also works on Windows, and LINUX despite its name). Much of what I do with MacMolPlt could also be done with Avogadro or Jmol (see below), but MacMolPlt is written specifically for GAMESS so it has some features lacking in the other two programs.

I usually do anything to do with structure in teaching with Jmol. Jmol is written in java, so the students can interact (rotate, zoom, download coordinates, etc.) with the models I post on the web. Furthermore, Jmol has a very powerful scripting language that allows one to create animations and control the animations via html. Note that virtually anytime you see a Jmol model on the web you have access to the coordinates by clicking on the window while holding down the ctrl key.

Finally, for anything to do with molecular dynamics I use Molecular Workbench in teaching. This wonderful program is molecular editor, MD engine, and analysis tool wrapped in one, and comes with a powerful scripting language. It is a bit like Avogadro/GAMESS/MacMolPlt/Jmol wrapped in one, but for MD instead of quantum chemistry. I'll show several examples of MW on this blog.

As you will see in future posts, I am a great fan of screencasts rather than written manuals and tutorials. Ironically, for screencasts I use a program that is neither free (it costs $100, ... ok $99) nor cross-platform (it works only on Macs running Leopard) called ScreenFlow. I am sure there are free alternatives, but based on what I read on google, this is by far the easiest to use. I can certainly testify that it is very easy to use.