Tunneling and STM

A few weeks ago, I gave a guest lecture (read: "I am at a conference that day, could you do it for me?") in a course entitled Unifying Concepts in Nanoscience. The topic was basic quantum mechanics (chapter 9 and a bit of 10 in Atkin's Physical Chemistry): particle in a box, etc.

These days, the first thing I do when preparing a lecture is to scour Molecular Workbench for useful animations, as I have discussed in a previous post. True to form MW did not disappoint, and I put together the following set of MW slides (note: you need to install MW first before clicking on it).

The screencast above shows how I used four of the slides to illustrate the concept of tunneling and and how it applies to STM.

Once again, I found animated simulations in general, and MW in particular, invaluable in bringing across complex concepts. And once again MW did all the hard work.

2012.09.01 Update: I made a few more screencasts of parts of my lecture

The Computational Chemistry Movie

Last week I presented the Computational Chemistry group to the new graduate students, and made the above movie for the occasion. As I mentioned in previous posts (here and here) I think molecular animation is a powerful but overlooked recruiting tool.

Most of the movie is one long Jmol script, but it took some post-editing of the resulting recording to smooth the transitions, i.e. remove the lag time when Jmol is computing the surfaces. The coordinates of the reaction and the large nanostructure is taken from Chemtube3D and this site, respectively. I have described how to extract the coordinates from a site using Jmol in a previous post. The orbitals and vibrational modes were computed using GAMESS and RHF/STO-3G.

The molecular dynamics animation towards the end is done with Molecular Workbench. You can find the simulation here, but you need to install Molecular Workbench to see it.

The equation animation and final credits are done with Powerpoint.

I hope you find the movie entertaining and useful, and feel free to use it (i.e. link to it, embed it, use in a Powerpoint presentation, etc.) if you do. Here is the page in blip.com where you can find the link and the source file in .mov and .flv formats.

Hearing voices

While looking for something completely different, I came across this online presentation by Martin Head-Gordon on computational quantum chemistry. The link is here and you click on View Presentation on the right.

The lecture is part of a series on modeling within nanoscience, and the hosting site, nanohub.org, has many other interesting things.

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?