Monday, June 15, 2009

The force is strong in this one

Here's a Jmol application I wrote to illustrate energy minimization/geometry optimization during a lecture in a molecular modeling course. Specifically, I wanted to illustrate how the atoms are displaced along the atomic force.

The arrows you see are the forces (or the gradient for the "start/gradient" button), which are the negative of the gradient (the derivative of the energy with respect to coordinates). The scale of the arrows relative to the scale of the molecule is arbitrarily adjusted to look nice. By clicking the buttons you are taken through the 8 steps it takes to reach a minimum. Select "wireframe" first to see the arrows more clearly.

Note that while the energy always goes down, the force vectors actually increase in the early steps. This is because if you simply displace the atoms along the arrows (i.e. perform the optimization in Cartesian coordinates) to decrease the angle, you also decrease the O-H bond length (which were fine to begin with). This increases the force on the atoms initially. However, as we get closer to the minimum the optimization fixes this.

The steps to create this Jmol application were as follows:

1. I performed the optimization with the GAMESS program at the PM3 level of theory (you can find the output file here).

2. I then wrote a python program (available here) to extract the geometries and gradients from the file and create the xyz+vib file that Jmol displays. The xyz+vib format is the normal xyz format (atom symbol, x, y, z coordinate) followed by the x, y and z component of the vector. You can access this file through Jmol as shown in a previous post.

3. I then created the html code, where I played around with scaling the vectors until everything looked nice. I finally decided on a scale factor of 10. You can access this html code as shown in a previous post.

PS: Very small arrows are not places correctly on the atoms. I think this is a bug in Jmol.


baoilleach said...

FYI, is a Python library that will extract the geometries from a GAMESS calculation. It doesn't currently do the forces though.

Jan said...

Thanks for reminding me about this. This might come in handy on something I am working on...

Geoff said...

You can display forces in Avogadro, although not (yet) from GAMESS calculations. It just requires some support in OpenBabel to treat the optimization as a conformer.

cclib is great, by the way.

Jan said...

Are there also plans to allow Avogadro to animate a geometry optimization file from GAMESS? I have used that in MacMolPlt and found it very useful for cases with poor convergence.

Btw I have submitted a bug report on the animation feature.

Geoff said...

Yes, there are definitely plans for that. :-)