Sunday, November 29, 2009

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.

Wednesday, November 25, 2009

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.

Saturday, November 21, 2009

Electron density

Fig3-1-1
Figure 3.1.
(a) A contour plot of the density of H2O 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 book and color figures

A big part of the motivation for this blog came from writing a book called Molecular Modeling Basics that will be published in May, 2010 by CRC Press. While writing the applications sections it was frustrating to turn the beautifully colored figures into black-and-white versions in order to keep the cost of the book reasonable. But it was also apparent that even colored figures in a book would be a somewhat poor substitute for the interactive versions they are based on. Especially, when turning them around to find just the right orientation for the figure. Wouldn't it be much better to have the reader decide for him/herself?

This is all a long winded way of explaining why there'll be a lot of posts with (color) figures that look like they came out of a book (they'll have figure captions below them). You can click on them for a bigger version. In many of the posts there'll also be a screencast showing how I made them, and an interactive Jmol version. They'll all be labelled "color figures from the book" so they should be easy find.

I always use a screen capture program, rather than saving a graphics file (see the screencast below for an example). I use a free Mac program called Snap and Drag, but I suspect there are plenty of other options for Windows and Linux.

Wednesday, November 11, 2009

A sense of scale


In chapter 1 of volume 1 of the legendary Feynman Lectures on Physics, Feynman starts by imagining zooming in on a drop of water, past amoeba and so forth, until one can see the water molecules. Starting this way is genius. The tiny length scales and the associated invisibility of atoms and molecules is the single largest barrier to developing a chemical intuition - a barrier that molecular animation can help overcome.

I am often thought of bringing Feynman's imaginary magnification to life by animation, so I as very happy when Nathan Baker brought this site to my attention. The above screencast shows it in action, but the real fun is interacting with it yourself. It's a brilliant piece of interactive animation.

The site brought to mind the granddaddy of them all, Powers of 10, which some kind soul put up on youtube.


Monday, November 9, 2009

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 Zn2+ 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.

Sunday, November 1, 2009

The trouble with transition metals I: molecule building



This post is long overdue. Quite some time ago (it has been so long I can't find the email anymore) someone asked me to do a post on transition metals.

Modeling of transition metal-containing compounds is notoriously difficult and will require at least two separate posts: one one building the molecules (that's this post) and one on SCF convergence problems.

The screencast shows how to build two molecules: a model of a zinc finger and a substituted ferrocene.

I build the zinc finger model in Avogadro and the only new trick is to pick the UFF force field, because that has parameters for bonds to transition metals.

However, I build ferrocene using Marvin Sketch and import the structure into Avogadro. This is because Avogadro doesn't have multi-center coordination bonds, while Marvin does. The "multi-center" points I add in Marvin Sketch show up as "dummy atoms" in Avogadro. As the name suggests these are not real atoms.

I then use Avogadro to add a functional group to one of the rings and optimize them while leaving the ferrocene part frozen. I could have added the group in Marvin Sketch but Avogadro gives me the opportunity to use autoopt to find the best geometry for the functional group. Note that because I am freezing the transition metal containing part I don't need to assign bonds between the Fe atom and the rings and therefore I can use the MMFF force field, as long as I delete the dummy atoms first.