Tuesday 28 August 2007

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.

13 comments:

Ψ*Ψ said...

I think phenanthrene should be less reactive than anthracene--2 aromatic sextets versus one. (Linear acenes have low HOMO-LUMO gaps in general.)
Very cool, though. You say it's possible to calculate these properties for pretty much anything? Hmmmm. I might play around with it. How long do these things usually take?

Egon Willighagen said...

Hi Felix, what license does the code have? It's important, so that people know if they may copy/distribute/change the code, or not.

Noel O'Boyle said...

Wow! :-)

The only other thing I can think to say is you might want to use something like glob.glob(os.path.join("..","*.mol2")) for identifying files matching a particular extension.

Noel O'Boyle said...

Actually, something else you might be interested in is to calculate the same HOMO-LUMO gap at a higher level of theory. You can easily do this in a few lines of Python if you install PyQuante (like the OpenBabel module, there's a Windows installer available).

Felix said...

thanks for the comments, nice to hear that it sounds kind of interesting. maybe that will make me work on it some more. that first thing will probably be graphics, like simple representations of the MOs

Ψ*Ψ: you can do calculations on π-systems with Hückel theory. with my package they can contain C,N,O (but adding more hetero atoms would not be very difficult).
the nice thing with Hückel theory is that you can get quite useful results very fast. the whole output for one molecule takes about one second. but you can still say something about reactivity when you look at the HOMO-LUMO-gap and regioselectivity when you look at the MOs (which are only numbers so far)

egon: I will add a license information when I upload a new. In principle: everything is ok but I would like to know about it

baoilleach: glob and join are nice commands (no program will ever crash because there is an error with slashes in the path)
PyQuante seems nice and I will probably take a closer look at it eventually. but I think I will stick with Hückel theory for this package. I guess if you want a higher level of theory you can skip my script and just use PyQuante

Noel O'Boyle said...

Regarding graphics, you may want to check out how PyMOlyze or ccp1gui do it.

Anonymous said...

Hello,
the code looks quite nice but I wasnt able to compile the python module on Mac OS X. With Linux I wasnt able to impot openbabel in python.
I will make a new try on the weekend.

I do not like python :(

Felix said...

all I know is that there was a problem in the OB-python-windwos distribution that is fixed now.
you may want to check the PYTHONPATH (if openbabel.py is not found) or BABEL_DATADIR (if the dlls aren't found) variables. or you ask the OB development team

Anonymous said...

I didn't fint the problem on my Mac, I think the bulid script for the python module is just broken but with Linux it runs fine now.
Looks great!

Ψ*Ψ said...

Just a quick comment--it would be much, much more useful if you added sulfur. ;)

Felix said...

anyway, nice to hear that you used it
There is no parameter for S in that list I used. They have halogens though ...

Desdicardo said...

Hi Felix,

Nice idea, with this code!
But the way, why it is no more currently available the code? When I want to see the code I cant find them at the link that is indicated in your post. I wander if it is possible to have a copy of this code from you because I want to play around with python and with this kind of application which involves science.
Thanks

Felix said...

Thanks. Now all the links should work again. I had to move the server and then I did not put the links correctly.

or directly to the download page:
http://homepage.univie.ac.at/~plassef6/chemprogs/python.htm

let me know if you have a question or if something does not work