Hückel package

Openbabel for Python is finally running on my computer. Here is what I did with it: a package that does Hückel theory for me. Thanks to numpy which does the math and pybel which reads the molecules, my package is only about 200 lines long (most of it comments and blank lines). I previously did something similar in Visual Basic. I had to do all the parsing myself. For matrix operations I took mathematica which I had to open seperately. And I could not distribute it easily. The bottom line is: Python is cool. This is a typical example of running the package with a mol-file of pyridine. The molecule is read in, information is extracted and the Hückel matrix is made. It is diagonalised and we get energies and MOs. The program also calculates the HOMO-LUMO gap and the mean π-electron energy if you want it to.

>>> import hueckel
>>> pyr = hueckel.hueckel('***/02_pyridine.mol', calc_everything=True, print_results=True)

Hueckel matrix
[[ 0.   -1.   -1.    0.    0.    0.  ]
[-1.    0.    0.   -1.    0.    0.  ]
[-1.    0.    0.    0.   -1.    0.  ]
[ 0.   -1.    0.   -0.83  0.   -1.  ]
[ 0.    0.   -1.    0.    0.   -1.  ]
[ 0.    0.    0.   -1.   -1.    0.  ]]

Energies
[-2.21192929774, -1.26984891306, -1.0, 0.748053194363, 1.0, 1.90372501644]

MOs
-2.21193 [ 1.10596  1.44632  1.       2.09318  1.10596  1.44632]
-1.26985 [ 0.63492 -0.19374  1.      -0.88095  0.63492 -0.19374]
-1.0 [ 1.  1.  0.  0. -1. -1.]
0.74805 [ 0.51933  1.      -1.38849 -1.26738  0.51933  1.     ]
1.0 [ 1. -1.  0.  0. -1.  1.]
1.90373 [-0.95186  0.81208  1.      -0.59412 -0.95186  0.81208]

pi electron number: 6

HOMO LUMO gap: 1.74805319436

mean pi electron energy: -1.35559273693

You can easily do calculations on a whole set of molecules. This prints the HOMO-LUMO gap and the mean π-electron energy for all the molecules in directory condensed_aromatics.

import os, hueckel
   os.chdir('***/condensed_aromatics/')
   file_list = os.listdir('.')
   file_list.sort()
   for file_name in file_list:
           if file_name.partition('.')[2] == 'mol':
               print '\n' + file_name
               h_mol = hueckel.hueckel(file_name, calc_everything=True, print_results=False)
               h_mol.print_results(mean_energy = True, HOMO_LUMO_gap = True)

In the output you notice that average resonance energy increases when you condense rings (this is a thermodynamical property I would say). But the HOMO-LUMO gap decreases making the substance more reactive. At 11 there is phenanthrene which has a much higher HOMO-LUMO gap then anthracene and should therefore be less reactive. I think I read that somewhere unless it was the other way around? And then there are a few more substances (with not really their IUPAC names) until we have a circle of 6 benzenes, no idea how that is called. The nice thing is that all is done in less than 5 seconds. It makes me think that numpy is cheating.

01_benzene.mol
HOMO LUMO gap: 2.0
mean pi electron energy: -1.33333333333

02_naphthalene.mol
HOMO LUMO gap: 1.2360679775
mean pi electron energy: -1.36832385059

03_anthracene.mol
HOMO LUMO gap: 0.828427124746
mean pi electron energy: -1.37955060707

04_tetracene.mol
HOMO LUMO gap: 0.589925798583
mean pi electron energy: -1.38504578655

05_pentacene.mol
HOMO LUMO gap: 0.439373742196
mean pi electron energy: -1.38836439167

06_hexacene.mol
HOMO LUMO gap: 0.338749064157
mean pi electron energy: -1.3906143081

07_heptacene.mol
HOMO LUMO gap: 0.268449344637
mean pi electron energy: -1.39225072443

08_octacene.mol
HOMO LUMO gap: 0.217562743235
mean pi electron energy: -1.3934980926

09_nonacene.mol
HOMO LUMO gap: 0.179638872553
mean pi electron energy: -1.39448163173

10_decacene.mol
HOMO LUMO gap: 0.150676414586
mean pi electron energy: -1.39527744692

11_phenanthrene.mol
HOMO LUMO gap: 1.09263880094
mean pi electron energy: -1.34260447994

12_benzophenanthrene.mol
HOMO LUMO gap: 1.06469427106
mean pi electron energy: -1.36802047201

13_naphtophenanthrene.mol
HOMO LUMO gap: 0.976434573677
mean pi electron energy: -1.38150272247

14_hexabenzo-circle.mol
HOMO LUMO gap: 1.06469427106
mean pi electron energy: -1.36802047201

If someone wants to try out that package, it can be found here. To use it you need Python, the numpy package and the openbabel for python distribution. And you can of course tell me if you like it. For everyone else: I can only recommend numpy and openbabel for python.

RMSD

RMSD (root mean square deviation) is the typical measure to compare different structures of a molecule. The nice thing for people with a linear algebra fetish is that the RMSD is nothing but the distance in the appropriate euclidian space. Let's see why, using this neat script for making LATEX formulae I found at A Zephyr in Time.

I don't know which indeces are least confusing but let's do it like this:
In structure A of our molecule with N atoms, the atoms have the following coordinates:



Structure B is the same molecule but bond lengths and angles are changed (or it is moved in space). The coordinates are



With the chosen coordinates the RMSD is defined as follows:



This reminds the attentive reader of the dot product she has been using since high school.



This is the length of the difference vector divided by , in other words the distance.

If you want some more linear algebra you can go on:

With two coordinate vectors a and b



we define the scalar product (using the Dirac notation to remind us that a typical overlap matrix element can be interpreted as a scalar product)



Now the RMSD between two structures A and B can be computed as the length of the difference vector between the two coordinate vectors a and b (containing the coordinates of all the atoms).



RMSD is brought down to something anyone can handle: a distance. And it is brought down to a clean mathematical construct, e.g. triangle inequalities immediately follow. I thought that was pretty cool.

Nice software?

I feel like blogging because it seems that my job is finally running the way it's supposed to. Let's review the software used (and complain).

There's the Python Macromolecular Library. It has tons of cool functions like telling you that your atom names are non-standard amino-acid-atom names or creating a chain with all your waters. It's a nice clean elaborate package but I think for doing so much work, they a left out important functions. For example changing fragment IDs is rather difficult or adding fragments to chains. And saving parts of a protein file in a new pdb is not convenient.

Babel and Open Babel are really nice programs but they are probably what drove me nuts most. All that changing between file types is pretty cool. But the problem is if you want specific output. Starting out with a simple parser would have probably been faster than all that messing with Open Babel and re-editing the files afterward. And where I really needed it, it failed me: for adding hydrogens. The interesting thing is that Open Babel protonated His and Arg, and Babel left both without protons. Both is of course wrong in physiological conditions. Using the pH-model function, which seems pretty nice, did not change anything.

Maybe the solution would be Open Babel in Python. You can probably do more when you are closer to the source. You can even Kekulize molecules there.