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.

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.

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.

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.

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.

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 Zn²⁺ 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.

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.