Saturday, July 25, 2009

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.

30 comments:

NUchem said...

Jan,

These screencasts are very helpful indeed. My computational skills are slowly returning to me after a decade of vacation :) Some more suggestions I've thought of could include..

1) TS calculations for an Sn2 reaction, or similar
2) NMR calculations

It appears that the literature is increasing with the number of NMR calcs for molecules. This is something I've never tried in GAMESS and did notice that there's an NMR option in MacMolPlt when generating the .inp file. Have you had experience with NMR calcs before?

NUchem

Jan said...

NUchem,

Screencasts involving TSs are planned in the not so distant future.

Unfortunately, the current implementation of chemical shift calculations in GAMESS is so inefficient that it is unusable for all practical purposes.

I have heard rumors that are more efficient implementation of NMR parameters is coming, but I don't know the time table. I'll try to find out.

The only free quantum software I know of that can compute NMR parameters efficiently is Dalton, but I don't have much direct experience with it.

I do know a Dalton input menu is planned for Avogadro, but I don't think they have gotten very far.

There is also ORCA, but I don't know how efficient the implementation is.

reve_etrange said...

Jan,
I've tried running these calculations on a couple of machines. The first time the energy returned after minimization in Avogadro was slightly high - imaginary frequencies were returned. Then, I perturbed the dimer manually during that step and got to a lower energy, leading to good frequencies. On my other system, the frequency calculation led to different numbers (and an energy that is a few kcals/mol higher). The electronic energy from the DFT calculation was the same in both cases.
I guess small differences in the starting geometry matter...but how do we know when our results are "correct?"

-da

PS Your site rules.

Jan said...

-da - If I read you correctly, you optimized the water dimer at the PM3 level, using 2 different starting geometries.

In the first case you eventually managed to reproduce the PM3 result of the post (i.e. a free energy correction of 13.82 kcal/mol).

In the second case you got imaginary frequencies, and an free energy correction that is a few kcal/mol higher.

However, if you compute the B3LYP/6-31G(d) single point energy using either geometry you get the same electronic energy (-152.7533 au).

Is that accurate? If so, I have no explanation.

Jan said...

PS - Thanks!

reve_etrange said...

That's what happened, except that when I got the imaginary frequencies I perturbed the system and got to real frequencies, with an energy that was still somewhat (~2 kcal/mol) than 13.8.

Does the geometry with imaginary frequencies represent the formation/breaking of the hydrogen bond?

-da

Jan said...

Imaginary frequencies for "optimized" geometries typically mean that the energy is not completely minimized with respect to some degree of free (the normal mode associated with the imaginary frequency).

Nilesh Tawari said...
This comment has been removed by the author.
Nilesh Tawari said...

Dear Dr. Jensen,
Thanks a lot for your useful blog and the more useful book. I am reading your book "Molecular Modeling Basics", one thing I can say that the book is simply excellent.
Thanks
-nilesh

Jan said...

I am very happy to hear that you find both the blog and the book useful. Thanks for letting me know; I appreciate it.

Anonymous said...

AMAZING man!!
Your blog is great!!

A question for you: I've performed some NMR calculations with GAMESS; can you suggest me a program that is able to open the output file, for visualize the NMR-spectra?

Thanks

Jan said...

Thank you!

I'm afraid I don't know know any program that can do this. Avogadro has a "Spectra" window, but I am not sure what it does. I suggest posting your question on a mailing list. Perhaps the Avogadro or GAMESS mailing list, or ccl.net

Anonymous said...

Suppose one was to spend a few days calculating something (a geometry optimization job) and then when gamess exited gracefully one was to press the "Clean Up" button in GamessQ prior to "saving output to folder".
How would one go about retrieving the output file in this case? (After finishing raving about how there should really be a time machine disk hooked up to this mac!)
"My friend" is having this problem is rather bothered because subgroup meeting is very soon.

Jan said...

Your files (sorry, your "friends" files) are probably in the spool directory: ~/.gamessqd. Check GAMESSQ preferences to be sure of the location.

Notice the . in front, which means the directory is invisible, so you have to use old fashion cd in the terminal window.

Anonymous said...

Hi prof Jan, I have a question: the energy value obtained after an optimization process using the pm3 method is in hartree or in kcal?

thanks for your reply

Jan said...

The energy marked "FINAL" is in Hartrees

Anonymous said...

unfortunately the .gamessqd was "cleaned up". Login in as root I found it where gamessQ pointed to. But there was nothing there :( except for the currently running jobs.
To all others who may share my anal retentive tendencies- be sure to save your output before cleaning up.. :)

Jan said...

Ahh, too bad. Thanks for reporting back on this.

Anonymous said...

hey. just starting out as a grad student doing some calc's on gamess. Love your site. after an optimization energy(elec) +energy(nuc.) is about -76 hartree. thats a hige number in kcal/mol.. what are the units displayed in the optimization output?thanks alot

Jan said...

The units are indeed Hartrees. It is the energy for the rxn electrons + nuclei -> molecule, which is why the magnitude of the energy is so large.

Paolo said...

Dear Dr Jensen,
let me tell you that thanks to your blog, I'm able to do useful things concerning my research job, with GAMESS obviously. I wrote in October asking about the computational approach regarding Diels-Alder reaction. Now my calculations are in a advanced state but I don't know why when I'm performing a calculation like an optimization using B3LYP/XXX, in the output file is not present the heat of formation. Instead, performing the same calculation with a semi-empirical method this value is present.
Now, is there a way to obtain the formation enthalpy using a B3LYP/XXX?

Thanks for the answer and for your availability!!

Jan said...

I am very happy to hear that you find the blog useful.

The semiempirical methods are parameterized against experimental heats of formation, which I why they print it out.

The heat of formation is computed based on the energies of each atom. This can be computed once and stored for the semiempirical methods.

But for ab initio methods the atomic energies would have to be recomputed for every combination of basis set and correlation method.

Paolo said...

Thanks for the answer!
And if I perform the energy using PM3/ with polarizations like (2d,p), is it legal? I tried that and the output file is valid, but I don't know if this way to do is elegant. I need to have enthalpy of formation to compare with energy calculated using B3LYP/6311G(2d,P).

Jan said...

"And if I perform the energy using PM3/ with polarizations like (2d,p), is it legal"

No. I suggest finding a molecular modeling book, and reading about PM3, to figure out why not.

You should compare methods by comparing their predictions (structures, relative energies, etc) of the item you are interested in.

Vicente said...

Dear Dr Jensen,

Congratulations for your good work done with this blogspot and your book. I needed information for installing Gamess and working with it, and the information in the blog and in the book has been very useful. I know that you have worked on PCM solvation model in Gamess and I would like to ask you a couple of questions related to that model. I have optimized the geometry of a Rh complex in water and ethanol with a model chemistry B3LYP/dzvp. I have used the line $PCM slvnt=water $end and the program has run fine. My problem is that I have to calculate solvation free energy of the complex,so I need cavitation and dispersion repulsion free energy. Adding Icav=1 I can do a single point energy calculation with optimized geometry obtaining cavitation energy, but dispersion repulsion energy are only calculated for Hartree-Fock not for DFT. My question is if disp-rep obtained from a HF calculation could be used for determining solvation free energy. Do you know a source of van der Waals radii for transition metals?. Thank you very much

Jan said...

I am very happy to hear you find both the book and blog useful.

The display-rep turn does not depend on the wave function type (as faras I know), so I am quite surprised it doesn't work for DFT. Anyway, you should be safe using the RHF value, as longs as you use the DFT geometry.

Wrt to vdW radii for TMs, I suggest you use ionic radii. TMs tend to have low electronegativity and are essentially ionized in most complexes. Ionic radii are easily found with google.

Manab said...

Hi Jan,

I am very new to GAMESS and would like to optimize my system of interest using RHF and 6-311++G(2d,p). Howeverm moment I use diffuse functions, I get error message.

Could you please let me know how can I run my GAMESS program.

Thanks and regards - Manabendra
(E-mail : param627@gmail.com

Jan said...

I suggest posting your question, with many more details, on the Google GAMESS group: http://groups.google.com/group/gamess?hl=en

my views said...

i want to know how can i calculate formation energy using gamess

deba said...

Hi Jan, I am very new in GAMESS(US). I want to know is there any special command to see the results of GAMESS output files..