Chemical Quantum Images

Excited states

2012-01-13T10:46:00.004+01:00

Our review article "Electronically excited states and photodynamics: a continuing challenge" just came out. Somehow I am even first author ...

The article is a summary about today's methods for computing excited states. And in particular there is a section about interchromophore interactions, which I wrote. If you are interested in any of that, give it a try ...

Orbitals in Jmol

2011-10-18T12:03:00.006+01:00

Actually automatic plotting of orbitals in Jmol is even easier than in VMD, considering that Jmol can directly read the output of many major quantum chemical programs, or alternatively Molden format files.

The only thing is that Jmol does not have all of its functionality easily available over a menu, so you should always have the scripting manual ready and use the console inside the program. But then it is rather easy to use. And the advantage of script based input is that you can naturally automatize repetitive tasks. For visualizing orbitals it can look something like this (you can just copy this directly into the Jmol console):


background white
mo fill
mo cutoff 0.03

mo homo

At this point you would probably play with the perspective a little bit. Then you continue by saving this picture and all other pictures of orbitals you are interested in, e.g.


write image png "homo.png"
mo homo-1
write image png "homo1.png"
mo homo-2
write image png "homo2.png"
mo lumo
write image png "lumo.png"
mo lumo+1
write image png "lumo1.png"

And you have what you wanted.

Actually the new Jmol lets you do even more and you can make linear combinations. For example if we combine the HOMO (14b_g) with the HOMO-1 (13a_u), we get an orbital, which is localized on one half of the molecule.


mo [1 183 1 182]

It is even possible to put up interactive models with Jmol. I used to do that, but then I lost my account on the server where I had this, and I did not want to set it up again. So nothing of that now ...

Orbitals, example

2011-10-12T14:44:00.006+01:00

Here is an example for automatic plotting of orbitals with Turbomole and VMD, as I explained in the last post.

What I am showing are the frontier orbitals of an adenine-thymine stack inside of DNA. The calculations were QM/MM calculations where the eletrons of the two bases were treated quantum mechanically and the remaining DNA and waters were considered as point charges. The orbitals are shown in order of decreasing energy with a qualitative label below them.

Thymine - π*

Adenine - π* (LUMO)

Adenine - π (HOMO)

Thymine - π

Adenine - π

Adenine - n

Thymine - π

Thymine - n

Automatic plotting of orbitals

2011-10-06T10:45:00.002+01:00

Back to computational chemsitry: I want to introduce a script that allows you to plot a number of molecular orbitals automatically in batch mode and yields nicely rendered image files without the need to sit there and wait for your program to compute every single orbital. All you need is VMD and a program that can make surface files (.plt, .cube, ...) for VMD to read. [1] In Turbomole, for example, you only have to add the line


$pointval mo  65-74

Then you would call the following script "vmd_get_plt.bash"


#!/bin/bash
# 1. call this script
# 2. open the molecular structure file in VMD
# 3. load the .plt files and some settings
#    - "Load state" load_all_plt.vmd
#    - click "Apply" in "Graphical Representations"
# 4. adjust perspective
# 5. "Load state" plot_all.vmd

out=load_all_plt.vmd
plot=plot_all.vmd
rm $out $plot

echo "axes location Off" >> $out
echo "display projection Orthographic" >> $out
echo "color Display Background white" >> $out
echo "menu graphics on" >> $out
echo "mol modstyle 0 0 Bonds 0.100000 10.000000" >> $out
echo "mol addrep 0" >> $out
echo "mol addrep 0" >> $out
echo "mol modmaterial 1 0 AOShiny" >> $out
echo "mol modmaterial 2 0 AOShiny" >> $out
echo "mol modstyle 1 0 Isosurface  0.035000 0 0 0 1 1" >> $out
echo "mol modstyle 2 0 Isosurface -0.035000 0 0 0 1 1" >> $out
echo "mol modcolor 1 0 ColorID 0" >> $out
echo "mol modcolor 2 0 ColorID 1" >> $out

N=0
for I in *plt
do
 echo "mol addfile $I" >> $out
 echo "mol modstyle 1 0 Isosurface  0.035000 $N 0 0 1 1" >> $plot
 echo "mol modstyle 2 0 Isosurface -0.035000 $N 0 0 1 1" >> $plot
 echo "render TachyonInternal $I.bmp" >> $plot
 N=$(($N+1))
done

This script will create two TCl scripts, which can be accessed by VMD through "Load State". The first one "load_all_plt.vmd" will load the plot files into VMD and set some parameters in a way that I think looks good. You can of course modify the initial bash script to change some of those parameters.


axes location Off
display projection Orthographic
color Display Background white
menu graphics on
mol modstyle 0 0 Bonds 0.100000 10.000000
mol addrep 0
mol addrep 0
mol modmaterial 1 0 AOShiny
mol modmaterial 2 0 AOShiny
mol modstyle 1 0 Isosurface  0.050000 0 0 0 1 1
mol modstyle 2 0 Isosurface -0.050000 0 0 0 1 1
mol modcolor 1 0 ColorID 0
mol modcolor 2 0 ColorID 1
mol addfile 11ag.plt
mol addfile 11b1u.plt
mol addfile 12ag.plt
mol addfile 13ag.plt
mol addfile 7au.plt
mol addfile 7b1g.plt

The second one is used for the actual plotting. First arrange the settings in a way that you like. Then call "plot_all.vmd" through "Load State" and it will create the image files.

mol modstyle 1 0 Isosurface 0.035000 0 0 0 1 1 mol modstyle 2 0 Isosurface -0.035000 0 0 0 1 1 render TachyonInternal 11ag.plt.bmp mol modstyle 1 0 Isosurface 0.035000 1 0 0 1 1 mol modstyle 2 0 Isosurface -0.035000 1 0 0 1 1 render TachyonInternal 11b1u.plt.bmp mol modstyle 1 0 Isosurface 0.035000 2 0 0 1 1 mol modstyle 2 0 Isosurface -0.035000 2 0 0 1 1 render TachyonInternal 12ag.plt.bmp

With knowledge of TCl programming it could of course be done in a more integrated fashion but I think it is convenient as it is. I will show some of the pictures soon.

Actually I found out that Jmol might be able to do it even faster because it can directly read the output files of quantum chemical programs and it computes the orbitals itself. And scripting in Jmol appears to be very straight forward. Maybe I will use that as well.

iCar

2011-09-18T13:21:00.004+01:00

Just a little thought for the day: Why are our cars not at least as smart as our phones even though lives depend on them?

It amazes me everytime how much smart phones can do. But why are cars as stupid as 100 years ago?

What's so difficult about Bluetooth communication rather than old fashioned horns and lights? Aside from being much more efficient it would also eliminate a large part of traffic related nuisance.

Why can a car not compute the ideal speed for reaching the next traffic light when it turns green? Why do people still hit the gas like maniacs in front of a red light just because it might turn green before they get there?

Is it so difficult for parking spaces to communicate with cars and tell them where to go? Do we really want all those cars doing extra loops with erratic speed changes and sudden stops, instead of a little bit extra wireless communication?

If you think about how smart phones are way beyond anything you would have imagined five years ago, you can understand how much creative potential there is also for cars. So why are we still crashing, honking and cursing like 100 years ago?

Waiting

2011-07-30T22:41:00.005+01:00

I have not had the motiviation to write too many science related posts recently. But maybe it is still ok if I post some arbitrary thoughts: It is interesting how much effort people spend only to wait longer.

Most prominently I notice that with flights: While I am enjoying the last minutes where I can sit comfortably and stretch my legs, I see how the queue in front of the gate is getting longer and longer. Apparently you are not going to get there faster if you board earlier. In typical airlines the seats are preassigned, which means that boarding earlier will not change anything. The only thing that changes is that you have more time waiting squeezed in the plane.

And leaving the plane is even worse: The seat belt sign is not even off when people get up, only to stand there bent over for the next 10 minutes. They are out of the plane a little bit earlier but what happens is that: we go into a last-in-first-out bus, which means that someone who just relaxed in their seat and let everyone else go first will be the first one in the airport. For some psychological reason it bothers me to see people do that and I feel like I should do that as well. But I should be happy that they let me get into the bus last ...

The other case is traffic lights. It is always striking me what kind of maneuvers car drivers tend to do only to get to a red light earlier, in order to wait there longer. On the one hand dangerous noisy lane merges on the large streets. And for me even worse is how car drivers sometimes pass bike riders right in front of red lights only to come to a full stop and cut off the bicyclist a few seconds later. And again this might even make you slower: if you are going slow enough that you do not have to stop at the red light, you will not only save gasoline. But it will be actually faster because you already have a non-zero speed at the time when the light is green which means that you have an advantage compared to someone who is standing there.

So either there is something really exciting about waiting that I have not understood yet. Or it is that our intuition sometimes gives us the wrong impulse in modern city life and we should spend more effort in active thinking.

Traffic

2011-04-25T02:58:00.001+01:00

"Traffic" by Tom Vanderbilt is a very interesting book that shows some of the not so obvious or intuitive facts about driving. I think his perspective is quite realistic and unbiased. And since as he said, driving is the "most dangerous thing most of us will ever do", it seems to make sense to take some time and think about some of those questions.

The main point is that using only intuition we will misjudge many of the risks and consequences of driving. One typical fallacy is that when people are asked how good their driving is, more than half of them will say they are above average. This is of course not possible,[1] people think that they are better than they are. A similar fact is that people in higher cars will tend to go faster. There is no logical reason for that. It is just a question of perception. A third point is that the design of the street has a significant impact on how we realize whether other people are affected by our driving.

The thing with driving is the high toll of people being seriously injured or even killed. Plane crashes and bus accidents are on the news, but car accidents are something where many of us have personal experience. Vanderbilt points out that even in high risk industries such a high rate of accidents could never be tolerated. Remember September 2001 - even in this very month more people died on US streets than from the terror attacks. I am always amazed by all the fire security measures - and I think we can really say that horror stories that you hear from earlier times are very unlikely to occur again. But why is our society not willing to reduce the amount of people injured or killed on streets? It is obviously possible to find a speed with no serious injuries.

Another point that he makes is how often car commercials tend to show dangerous driving. That is almost like a beer commercial would feature a drunk driver. From my analysis of US car commercials I can say that some other themes might be intended to apply to the herding or pecking order insticts. But practicability or comfort seems to have only a minor role. And even more usually it is assumed that the person watching the commercial will take up a loan over five or six years just to pay this car - like that person does not have anything better to do with their money than paying for a car which is a little bit bigger and more noisy than usual. I guess the reason is prestige. But it is of course sad that society gives up so much for prestige. Isn't there a way of showing the status in the pecking order that does not annoy anyone else?

With respect to smoking there was a lot of increase in awareness recently and people are realizing that smoking may seriously affect other people not smoking. It is not perfect yet, but at least I am noticing that smokers would sometimes apologize for smoking in the presence of non-smokers. But I don't see that with cars. I always see cars racing across pedestrian crossings as soon as it looks like the person waiting is scared enough to give up his right to cross. People have realized that smoking is dangerous and a person smoking a really strong and filterless cigarette will seem like a nut rather than a really cool person. But for some reason someone driving an unreasonably large car will actually have a high social status, even though he is scaring little children off the street, wastes precious recources, and pollutes everyone with noise. I think here society still has to wake up.

[1] Precisely, one would have to specify that the median is meant, not the mean. If the distribution had a high skewness, i.e. there were some really bad drivers, then it would in fact be possible that more than half are above the mean.

DNA base pairing

2011-04-04T03:17:00.010+01:00

Another very interesting point is the geometry of Watson-Crick base pairs in DNA.

First a look at the adenine-thymine base pair. Thymine is on the strand whose 5' phosphate group is shown on top. The strand with adenine goes the opposite way and you see an open 3' oxygen. There are two hydrogen bonds in the AT pair.

Let's look at this base pair from the other side. The important point is that the backbone looks the same as in the above image. And ths is why you can have alternating A's and T's in a strand. Apparently this is what Francis Crick pointed out to James Watson, when Watson had suggested base pairing. Considering that this had to be pointed out even to the person winning the Nobel prize for discovering the DNA structure, I think it is safe to say that people without this understanding of helical symmetry do not have to feel bad about it. But apparently it is a non-trivial fact.

The next non-trivial fact is that the backbone of the guanine-cytosine pair looks just like the AT and TA pairs. Here you can see the base pair with its three hydrogen bonds and which is thereore bonded a little bit more strongly.

And again, CG also fits.

It seems therefore that one reason for the specific choice of bases is that they all fit into the same helical structure. And therefore it is possible to have a stable helix which contains the genetic code through alternating bases. Another source of selection pressure was probably photostability, as explained here.

DNA close up

2011-03-29T04:00:00.010+01:00

I wanted to look a little bit more at the details of the DNA structure. In part motivated by James Watson's book (described here). Some time ago I showed pictures of the helix, but now I want to show some molecular details.[1]

This is the backbone of DNA in its B-form. The phosphate is in front. The 2-deoxyribose is the vertical 5-ring. Its anomeric center is bonded to the nitrogen of the base and the oxygen closing the furanose ring. The 3' oxygen bonds to the next phosphate group. The 5' oxygen leading to the other side is bonded to the carbon outside of the ring.

In principle a pentose has five active groups. In DNA all are saturated and no OH groups remain. Why?
1' bonded to the nitrogen
2' that's the "deoxy" in "deoxy-ribonucleic acid"
3' bonded to phosphate
4' forming the ring
5' bonded to phosphate

What happens if we put an OH group at the 2' position?

The new OH groups is shown with its full van-der-Waals radius. There is clearly not enough room for the extra oxygen atom, considering for example that the other oxygen atoms should be just as big. It is in fact not possible to push the atoms out of the way in any simple manner. Therefore RNA does not form a B-helix.

The structures were created with the nucleic functionality of Tinker. The pictures were drawn with pymol.

[1] One thing I really learned in science is quoting myself ...

Traveling to Texas

2011-03-06T01:03:00.002+01:00

It is my third time at a research institution in the USA, this time at Texas Tech. The first good impression was that when I filled out my "Electronic System for Travel Authorization" form, they told me I did not even have to fill it out because I was already authorized - and that even though last time I went to the US there was not even such an electronic system. That's good administration. When I got there this somewhat changed after waiting in line for an hour while my connecting flight was leaving. Even worse was that when it was finally my turn I was a little bit detained because I did not know if I was "visiting" or "researching". Why would you do all the J1 bureaucracy if you are just visiting a group to make some personal contacts? Anyway the second guy told me that it was no problem and I could enter.

What I really like about US universities is that they are trying to make them nice and productive places. The campus looks nice and when you go to your office it looks like it is all prepared with the idea in mind to make an open environment for good research. Back in Vienna we are in a building that was planned to be torn down 20 years ago. For example we do not have enough places to put up our posters, which I always think is sad. The only thing they did to improve our work is putting in new fire doors and some ugly new pipes for fire extinguishers. I am not saying it is not necessary to have protection against fires but I don't like the idea that everyting is built upon formal fire protection rules rather and no-one even thinks about what kind of measures would be important to provide a productive research environment.

The US mentality toward research is something that really motivates me to come here. But then some other things kind of strike me. Here it is mostly the huge abundance of SUVs and pick-up trucks. It makes me wonder what it is that people are trying to compensate ... But it's ok I guess, if for some people a stretch limo is the only way they can find a mate, aside from dangerous penis enlargement surgery, then I guess it is socially just that stretch limos exist. Or is it more the mentality of the type: "With all those crazy people out there driving SUVs, you'd better get the biggest car that you can!" Anyway, on most of the side streets car drivers are much more polite than I am used to in Vienna, so it's ok.

The step that follows after traveling is of course reimbursement, as described nicely at PhD comics (1,2). You give a free loan, collect tons of paperwork and then you carry it to the financial department on your knees. We even have to show our credit card statements - my opinion is: if they want a credit card statement, they should also give me a credit card and everything is solved. But what can I do. Interestingly the hosting US institutions were always much more generous, insisting on paying some extra food money, without really so much bureaucracy. But I think guests at our department at home also do not have to do all the paperwork. So it seems like departments always treat their guests better than their own people.

Drugs in our food

2011-02-10T17:15:00.010+01:00

"Opium fürs Volk" ("Opium for the people") by Udo Pollmer was a really interesting book showing how many different drugs and psychogenic substances are in our food. There are of course a lot of substances in plants which the plant needs to be healthy and which are also healthy for us, like vitamins and nutrients. But there is a whole other set of substances which plants or mushrooms use to fight enemies, or substances which incidentally interact with molecules in our body. And this second category is what this book is about.

Everyone is aware of things that are clearly toxic, like deadly nightshade or a toadstool. But the point is that as soon as things are not so toxic anymore we start to like them - as drugs. One prime example is the toadstool "fly agaric". It is most probably what sent Alice to Wonderland. It may also be what caused the ecstasy in the Dionysus cult. It was probably a major ingredient in Asterix' magic potion, causing the frenzy "Berserkerwut". And for some reason in Austria we give each other little red mushrooms with white spots for the New Year, to bring you luck - next to the more common four leaved clover leaves.

Another question this book helps to answer is what is so special about chocolate. Why doesn't just mixing fat and sugar do it? The answer is at least seven different kinds of substances, coming from cocoa, the production process and being built in our body. Here I want to focus on salsolinol.

A major effect of this molecule is as a monoamine oxidase inhibitor. This means that monoamine neurotransmitters, e.g. serotonin and dopamin, are not broken down. But they are transported to our brain and make us feel good. Salsolinol is also what explains the popularity of bananas. It is found in particular in the brown spots - maybe that is the reason why banana cake is usually made from really brown bananas. Aside from the psychotropic substances in hops, salsolinol (formed in the body) may even be another reason for the popularity of beer.

Maybe I will write about some more things in another post. Cola, why is it so much more popular than the other drinks containing sugar, caffeine and artificial coloring? Tomatos - "Paradeiser" (paradise fruit) in Austria and apparently also in the Czech republic, "Pomodoro" (gold apple) in Italian. Spices - why are the so much wanted, in particular Safran. Food processing - what is the deal with Aged 18 Whiskey, Balsamico Vinegar, barbecueing, ...

Charge transfer

2011-01-24T13:00:00.011+01:00

My new paper just became available online (DOI: 10.1063/1.3526697). The purpose was to present and test an ab-initio approach to simulate charge transfer dynamics. Aside from that we discuss some of the theory behind it and how it applies to direct non-adiabatic simulations, in particular the relation between Landau-Zener theory and Marcus theory. As a model we looked at hole transfer between two ethylene molecules bridged by up to three formaldehyde molecules.

Of course the next step is to extend the approach to more interesting systems. But the nice thing was that most of the physcially important parameters could be adjusted easily. So it was a good way to learn the basic phenomena of charge transfer.

I showed in a previous post what such a dynamics looks like. Here I want to talk a little bit about some of the math behind it. Before we did this work, I never really understood what a non-adiabatic coupling vector was. I will briefly describe that now - for a longer explanation you can of course download and read our paper.

The components of the non-adiabatic coupling vector
are obtained as the derivative of the electronic wave function of one state projected onto the wave function of another state. R_j means the displacement of a nuclear coordinate. Note that this is a derivative in the Hilbert space of electronic wave functions, parametrically dependent on nuclear coordinates.[1]

The interpretation comes when using the mixing angle η between diabatic and adiabatic functions. Because it turns out that
What I am doing now, is to start with a delocalized charge (η=π/4) and localize the charge (η=0) by changing the bond length alternation.

At an intermolecular distance of 5 Angstrom, the interactionbetween the fragments is quite large. The wavefunctions change slowly over a wide geometric range and you have a wide peak of the coupling vector accordingly.

At an intermolecular distance of 7 Angstrom, the interaction is almost vanishing. The wave function localizes after only a small geometric displacement. Accordingly there is a highly peaked coupling vector.

According to the second equation above, the area under both curves is π/4. If you look at the graphs you can see that this could be right, i.e. the area should be about unity.

[1] Similar difficulties occur in the derivation of the Hellman-Feynman theorem - in that post there is even the third component λ which represents the possibility of having a not converged wave function...

No slides, no sympathy

2011-01-13T21:24:00.008+01:00

Is it rude to walk out of a boring talk or a social responsibility?

Many things can happen with a talk: the accoustics may be bad, the speaker may be nervous. It's even ok if he/she is not so firm with every underlying detail. But for me the borderline comes when the speaker does not even care enough to prepare slides.

I was really excited about a public talk about Schrödinger "50 Years After". I just like the discoveries in those times and the stories behind them. And it seemed to be a possibility to show my girlfriend a little bit of what I am doing. I even cut track practice short to get there on time. But what happens: This man just starts reading in a monotonous voice. People start falling asleep. No slides or anything exciting to wake them up.

It is a petty in my opinion: It would take about five minutes to prepare slides that contain the names of the scientists mentioned. 10 more minutes to get their photos out of google. Maybe an hour to prepare some graphics. Then people could at least take a little bit home. And he does not even do it himself - a lot of people (including me) would be happy to see things they prepared in a big auditorium.

So I just walked out when I noticed that my girlfriend was falling asleep. It is not so much that I minded sitting there for another half an hour. I just think it requires a statement if a speaker does not even care enough to prepare slides.

Atkins

2011-01-06T19:22:00.007+01:00

The book I am currently reading is Robert Atkins' "New Diet Revolution". It is not that I am trying to lose weight but I just wondered what is behind it. And I think he has a lot of good points that standard nutritional theory has not fully appreciated yet.

The main claim of Atkins is that sugar and refined carbohydrates are more problematic than fat. This goes for both health and weight issues. I did not realize this initially but it makes sense. To make an analogy: flooding the organism with sugar is like pouring oil into a fire - you are adding highly reactive fuel. In more chemical terms: glucose being an aldehyde is a reactive molecule. And if too much is present it will cause problems, for example react with the amino groups of proteins - the Maillard reaction. What does the body do against it? Produce insulin to make sure all the energy is stored in a safe place - as glycogen or fat. And if everything is stored, there is probably enough excess insulin to lower our blood sugar level enough for the next hunger attack and cravings of carbohydrates. Aside from increasing fat tissue, too much insulin production may eventually also overstrain the pancreas and cause diabetes.

Another interesting point is that food energy, i.e. "Calories", is not everything. The idea is that the standard oxidation energy determined in a bomb calorimeter cannot necessarily be transferred 1:1 to energy available for the body. Atkins states that he has examples where a person would gain weight when consuming a specified amount of calories in carbohydrates and lose weight when they would consume the same amount of calories in fat. The energy conservation law also holds in biology of course, so where does the energy go? Probably into body heat. This leads to what another interesting nutritionist, Udo Pollmer, says: That people on a diet (which is probably low fat) feel cold quicker. [1]

As I just discussed, calorie counting (which seems a natural application of the energy conservation law) may not even make sense from a purely biochemical point of view. A point which is probably more important, that both Atkins and Pollmer mention, is psychology. Basing a diet on severely restricting yourself on such an abstract thing as the number of calories does not work for most people. The body has a built in mechanism to control food balance and it seems difficult to fight that with willpower. Eventually one would just stop the diet and "binge eat". Or maybe your body will get the calories it wants in a more subtle way, for example through high calory drinks (cola, fruit juices, alcoholic beverages, milk, ...). So what restricting calories does, is probably just telling your body that food is scarce and it should try to get as much as possible from wherever it can.

[1] edit: actually another reason is that in ketosis energy is lost in chemical form, i.e. high energy molecules like ketone bodies are excreted.

The Double Helix

2011-01-04T17:02:00.006+01:00

A very interesting book I just read is James Watson's "The Double Helix", his account of the discovery of the structure of DNA. It is a very fun to read story, showing how they finally came up with the double helix structure and the famous Watson-Crick base pairing. It also gives a lot of examples of how science should or should not be done. It is intersting to see how inefficiencies in the system and personal hostilities had to be overcome before the discovery was possible. He has a very interesting personal and funny writing style but is in some cases a little bit harsh, especially against co-researcher Rosalind Franklin (for which he apologizes in the epilogue, though).

I think the key to Watson's success was interdiscplinarity.[1] Only by bringing together x-ray data, experiments from organic chemists, and his own model building was it possible to solve such a complex structure with the intstruments available in the 1950s. Francis Crick is described as a genius, but as a character who most people could not deal with easily. Rosalind Franklin appears as an excellent cristallographer but threatened by a male dominated world. Sir Lawrence Bragg was the head of the institute who could see how his famous equation lead through this breakthrough discovery.

The race starts with Linus Pauling's discovery of the protein α-helix. After that, Watson and Crick as well as Pauling come up with a triple helix for DNA, with the backbone inside. Franklin has some clear evidence that the backbone should be outside but does not open up her data for display. Eventually Watson tries out models with the bases inside and the backbone outside. It is stimulating to hear of his excitement as he develops the concept of base pairing. For his success he had to collect information from quite a number of people. Franklin and her boss Maurice Wilkins provided the x-ray data. Crick was able to derive and solve the diffraction equations for helices. The work of Erwin Chargaff provided the vital clue that the amounts of adenine and thymine, as well as, cytosine and guanine, were equal. A direct connection to structural chemists was important to get the correct tautomers of the bases (which were wrong in many textbooks at the time). Watson's final discovery was that he could form AT and GC base pairs which had about the same shape - therefore yielding a regular helix. To compare the results to the x-ray data they had to create a molecular model. At that time that meant putting it all together into a large metal model after asking the workshop to make the pieces. Finally everything worked out and they could publish their paper with the nicely understated final sentence:

It has not escaped our notice that the specific pairing we have postulated immediately suggests a possible copying mechanism for the genetic material.

A point worth mentioning is the secrecy. It seems that the structure could have been found much quicker, had Wilkins and Franklin released their data earlier. I don't know if this is different in today's "publish or perish" situation, that many results are released as soon as they come. And that they are not kept in secret for so long. But maybe for the big things it is still the same.

Finally it was interesting to read how Watson seemed to be living a very easy going life while he made the most crucial discovery of base pairing. Maybe I should tell that to my boss more often: "When Watson invented base pairing, he played tennis every afternoon. Bye!" At least as far as creative break through work is concerned, flexibility seems to be very important. That is what Lee Smolin points out in his book "The trouble with physics". But I guess a large part of my work is just related to more or less standard data production. That could mean it would be better if I go to work at /:30 and wear a tie...

[1] Interdisciplinarity is a little bit overly popularized these days but I think it is still an important feature to be able to look beyond the narrow horizon of one specific field.

using "tee" for interactive input generation

2010-12-17T14:27:00.008+01:00

Here is just a little rather simple thing which took me a long time to find out.

The question was: If you have a program that needs interactive input generation, can you still create something like an input file? The motivation was mainly the define program of Turbomole. Assuming you have 100 molecules (or structures of one molecule) and you want to do the same type of calculation on all of them. Running define for 100 times is not really the most exciting option. So what can you do?

You run
tee define.in | define
on the first molecule and use the interactive input facility. The file define.in will contain your input.

For the next 99 molecules you can just use
define < define.in > define.out
in any kind of loop structure (using separate directories for each job).

In that fashion you had the advantages of an interactive input generator but you are still able to automatize things.

DNA photostability

2010-12-01T09:53:00.009+01:00

It's been a long time since I have written anything. Too many other things to do. But now I had the chance to play around with graphics, work related. This graphic shows the response of UV bases to sun light. It was for a promotion (DE, EN) of this PNAS article, which explains the relaxation mechanisms of DNA bases after UV irradiation. I do not have anything to do with the paper but it is still nice to see my image on the front page of the homepage of the University of Vienna.

The images correspond to minimum energy conical intersection structures of the four DNA bases. At such geometries the ground state energy is strongly raised but the energy of the biradical excited state is not so much affected or even lowered. Therefore ground and excited states can become degenerate and a transition can occur. It is probably through such geometric distortions that DNA bases can efficiently relax to the ground state on a pico second time scale. Then only a slight chance of undergoing any photochemical modifications remains.

Through this inherent photostability in can be assured that only a small fraction of DNA bases undergoes photodamage. The small part that does react, for example to pyrimidine dimers, can be repaired by enzymes.

Now the technical part about the image: The background comes from a picture of the sun against a blue sky in the Austrian alps. The molecular structures were rendered with pymol. The trick was to use a transparent background in the rendering

set ray_opaque_background=0

Through adjusting the light setting I tried to get the light source to where the sun would be for the separate structures. e.g.

set light, [-100,-20,0]

to move the light source to the left and a little bit lower which would have been the appropriate setting for the adenine on top.

By reducing the ambient setting I tried to make this more pronounced. For this I used the Setting/Edit All in the menu, which is in some cases more convenient than the command line.

Recursivity

2010-04-20T19:16:00.004+01:00

Since I just mentioned the topic: here is a nice webpage illustrating the concept of recursivity.

(The only problem is that the author is a little bit too excited about his stupid joke.)

Self Rep

2010-04-10T19:47:00.002+01:00

What is the output of this python program self_rep.py?


def self_rep(a): print a; print 'self_rep('+chr(34)+a+chr(34)+')'
self_rep("def self_rep(a): print a; print 'self_rep('+chr(34)+a+chr(34)+')'")

Try "cat self_rep.py", "cat self_rep.py | python", "cat self_rep.py | python | python", ...

This self repeating program is a nice idea I found in the book "Gödel, Escher, Bach" by D.R. Hofstadter. In general the book deals with how Gödel's incompleteness theorem Escher's drawings and Bach's music are related to the notion conciousness. By self-reference and self-repetition.

Gödel's string in plain language can be formulated as: "This sentence cannot be proven in number theory". It is apparently true (if it would be false it would be possible to prove it which means it would be true, which is a contradiction). Hence there is a true statement which cannot be proven. The other question is how this string can be a part of number theory - for that Gödel introduced his Gödel numbers. Finally, with a construction like the program above, he could have the string point at itself.

To finish this post, here is another program:


def self_inc(i,a): print a; print 'self_inc('+str(i+1)+','+chr(34)+a+chr(34)+')'
self_inc(1,"def self_inc(i,a): print a; print 'self_inc('+str(i+1)+','+chr(34)+a+chr(34)+')'")

Now you can try it again: "cat self_inc.py", "cat self_inc.py | python", "cat self_inc.py | python | python", ...

Charge transfer dynamics

2010-03-16T17:46:00.005+01:00

It's about time for me to put up another post here. What I chose as a topic is to present little bit of the work that I have been doing over the last year. Of course it went a lot slower than we thought it would but at least we already understand what the problems are...

And I can already make movies that I think are pretty cool. What I am showing here is the charge transfer dynamics of the ethylene-formaldehyde-ethylene complex radical cation. This is of course a somewhat artificial system. But it is good for understanding the physics and the computational problems related to such processes.

The main point of interest is the non-adiabatic coupling between electronic and nuclear degrees of freedom. The Born-Oppenheimer approximation breaks down because electron tunneling between the ethylene molecules happens on the same time scale as the nuclear dynamics. This phenomenon is traditionally described by Marcus theory. Here we are applying semi-classical surface hopping dynamics, as provided by the Newton-X program package.

In the movie the charge is shown in red. It can be seen how it hops between the ethylene molecules. But not every passage of the transition region leads to a charge transfer. It may be inhibited by non-adiabatic coupling. In these simulations this is represented through a "surface hoppping" to the excited state. In the video it is shown as a yellow background.

For comparison I am also showing the static picture here. Plot (b) shows the location of the charge and can be compared with the red color in the above movie. (a) is the energy gap and it is marked if the trajectory is in the excited state. (c) shows the CC distances of the ethylene molecules.

TwirlyMol (2)

2009-11-27T20:00:00.010+01:00

It looks like I cannot show more than two molecules on one page. Here is the third one I wanted to show (you may have to click the header of this post to see it):

For molecules that are not in their data base I can enter the molecule as InChI (which I can create through open babel).

This is the molecule of my Master's thesis [2,2'-bipyridyl]-3,3'-diol with the InChI 'InChI=1/C10H8N2O2/c13-7-3-1-5-11-9(7)10-8(14)4-2-6-12-10/h1-6,13-14H'.

[To correctly translate it into a form that the browser interprets correctly, I can convert it on their homepage (I enter the InChI and take the converted one from the link it produces). Then it looks like this: 'InChI%3D1/C10H8N2O2/c13-7-3-1-5-11-9%287%2910-8%2814%294-2-6-12-10/h1-6%2C13-14H']
[Edit:] The problem was that I did not put the 'InChI=' prefix and it is not necessary to do this conversion.

In fact the molecule should be all in one plane and have intramolecular hydrogen bonds but it is still cool that the tool can create the structure from just one line.

TwirlyMol

2009-11-27T18:14:00.014+01:00

The Chemical Identifier Resolver-TwirlyMol tool (for more information: Noel's post) is pretty cool. All I have to do is enter two lines of code to get an interactive look at my molecule.

For example this would be 'adenine' (the tool recognizes the word and produces the correct structure)

Or 'acetylcholine':

So far I have been using jmol for these things. jmol is more powerful but it also takes much longer to load and it has some problems when several applets are open on the same page. And for me there is the problem right now that all my applets will stop working when they delete my account at my old university (and I don't know if I feel like setting it up somewhere else and rewriting the links). Check them out a last time at the images label. On a bad day firefox may freeze if you do so ... but it is still worth it.

Density (2)

2009-11-16T15:59:00.005+01:00

For a long time I did not understand the connection between density matrices as shown in my other post and density matrices as they are used in quantum chemical programs. On first sight it is not quite apparent that you can get the matrix representation of the 1-particle density matrix in an orbital basis in the following way:

The operators in the second expression are the creation and annihilation operators from second quantization. The expression means: you take an electron out of orbital j, put it in orbital i, and overlap this wave function with the original wave function. Note that the first scalar product is in 1-particle space, the second one in n-particle space.

To reach this conclusion you just have to write it out explicitely and rotate the integrals a little bit. With the definition from the other post, it looks like this.

As a next step you also write out the n-1 integrations that are involved in the density matrix definition.

You rearrange it a little bit:

The crucial point is the interpretation of the terms in the parentheses. To do this, we expand our wave function in Slater determinants constructed from an orthonormal orbital set including φ_i and φ_j. Formally we can construct any wave function in this way in a complete CI expansion.

The question reduces to: How does multiplying with an orbital and integration change one such n-particle determinant? - If the determinant contains the orbital, one obtains the (n-1)-particle determinant where this orbital is taken out. If it does not contain the orbital (and therefore only orthonormal orbitals) one obtains 0. This operator is usually called annihilation operator and written as a_i.

We may rewrite this equation with the adjoint operator.[1]

This offers an expression that can be efficiently evaluated. In a determinant basis all you would have to do is check whether the orbitals are present for every pair of determinants - no explicit integrations or things of that sort are needed.

[1] The adjoint operator is usually called creation operator. It adds an orbital to a determinant or makes it 0 if the orbital is already there. The explicit definition would be something like multiplying with the orbital and antisymmetrizing the expression.

Another interesting point here is that we have scalar products (= matrix elements) with different numbers of particles. For every k, the k-particle space is constructed as a Hilbert space, i.e. you can form scalar products. It is not possible to have a scalar product between wavefunctions of different particle numbers. Therefore the direct sum of all k-particle spaces (called Fock space) is indeed a vector space but not a Hilbert space. edit: you can construct the Fock space as a hilbert space if scalar products with different numbers of particles are just zero.

Hellmann-Feynman theorem

2009-10-26T14:27:00.010+01:00

If the Hamiltonian of a system changes, the change in energy should be related to both the change in the Hamiltonian and the response of the wave function. The Hellmann-Feynman theorem states that for an exact eigenfunction, only the derivative of the Hamiltonian is decisive. This statement can be generalized to any kind of variational wave function. And since I finally kind of understand the math behind it, I will put up some formulas. The derivation is similar to this Ref. [1].

We start with some parameters describing our wave function, (e.g. orbital coefficients and expansion coefficients ).

And we need some kind of a rule to make an n-electron wave function from λ at a geometry R [2].

We can compute the energy of this wave function in the usual way

We may formally view this as a function of λ (not any more a bilinear function) and rewrite it according to

This way all geometry dependence is in the Hamiltonian - the actual geometry dependence that comes from the motion of the nuclear charges and the geometry dependence that arises from the wave function expansion (e.g. motion of atom centered basis sets).

The above expression gives you the energy for any kind of expansion λ at any geometry R. Now in practice we will have some kind of a rule (or mapping) on how to create the wave function (e.g. MCSCF) at every geometry.

Now, the geometry dependent energy with respect to this rule is given according to:

Next we can compute the derivative of this expression with respect to a coordinate x in R (derivatives will be denoted with a superscript (x) or (λ)). This requires the use of the chain rule [3]

or with the definition from above

The first term evaluates the change in the Hamiltonian with the reference wave function. The second term is the gradient of the energy with respect to the expansion parameters multiplied with the change in expansion parameters.

If λ was determined in a completely variational way, then

The corresponding second term cancels and we arrive at the Hellmann-Feynman theorem.

The significance is that with a completely variational method, we do not have to worry about the response of the wave function to the geometrical perturbation and the formulas simplify. In real life MCSCF (and the special case SCF) is the expansion where this holds. For such a wave function the gradient comes directly from combining the density matrix (from last post) with the atomic derivative integrals. CI gradients are more complex because the orbitals are not optimal with respect to the CI expansion coefficients. And in coupled cluster not even the amplitudes are variational.

[1] R. Shepard. The Analytic Gradient Method for Configuration Interaction Wave Functions. in Modern electronic structure theory. 1995 ed. D.R. Yarkony

[2] Note that there are three types of coordinates: wave function coordinates λ, coordinates of the Hamiltonian R (I am thinking here mostly about the nuclear geometry but it could be anything else as well) and implicit electronic coordinates that I am not writing out.
- The λ in the wikipedia article corresponds to my R actually (I hope it is still understandable).
- A slightly tricky thing in this context is that in general the same λ can give a different wave function at a different geometry. This happens because atom centered basis sets move around in space with the atoms.

[3] The notation is kind of difficult with all these different types of coordinates (maybe not written in the optimal way here) but in principle this is a straight forward application of the chain rule. You could also do it with Taylor expansions truncated to first order - but chain rule is prettier ...

Density

2009-10-20T19:00:00.001+01:00

Density matrices are not only nice for mathematical reasons, they also offer a stringent physical definition of orbitals. The point of this post is to show how you can create what is called "natural orbitals" from any kind of wave function. These things are also interesting with respect to orbital imaging (as introduced here and discussed for example here).

In order not to make any prior assumptions, I have to use general wave functions. In real life you would usually look at these things from the perspective of an orbital basis. And then the expressions would be simpler (maybe I will show that, too).

The one particle density matrix γ₁ of an n-particle wave function Ψ is defined as:

This quantity is a "convenient partial sum". It lets you evaluate all the properties of 1-particle operators without the need of the whole wave function. This operation (in a finite orbital basis) is carried out in many quantum chemical codes.

Because the particles are interchangable you only need one such density matrix. (An n-particle density matrix which is analogously defined, lets you evaluate the matrix element of any n-particle operator.)

First we have to define what a one particle operator is. I think the most straight forward (and hopefully correct) definition ist the following:
An operator A is said to act only on r₁ if

The idea is that you can do all the integrations with respect to the other variables ahead. You need some variable renaming to make sure that the operator only acts on the second function.

This contains already the density matrix

From the density matrix we can construct, what is called "natural orbitals" without the need of any a priory definition of an orbital basis.

Maybe it would be easier to understand it in the matrix representation but in the general case we can also define an operator of the following form.

We may also symmetrize it (again easier in cartesian coordinates; edit: the 1-particle density of a wave function is already symmetric). This symmetrized operator has clearly defined orthogonal 1-particle eigenfunctions, the so called "natural orbitals" (and occupation numbers as eigenvalues). The natural orbitals represent the 1-particle density matrix. And this one particle density matrix contains all the information needed for 1-particle operators (like e.g. the dipole). So natural orbitals are something kind of physcial in my opinion.

The eigenvalues of the density matrix are the occupations. In a Hartree-Fock wave function the natural orbitals correspond to the usual HF orbitals and the eigenvalue of occupied orbitals is one (zero for virtuals). The density matrix is idempotent (or equivalently a projection). In the general case there is no clear distinction between occupied and virtual orbitals and there is strictly spoken no such thing as a HOMO or LUMO. But then we notice that HF represents most wave functions quite well in "well behaved" cases. So you may still think about HOMOs and LUMOs even in a more strict sense.

cgi

2009-10-15T18:35:00.005+01:00

This page is my contribution to a world without spam.

The idea is that it gives you a new random email address everytime you visit it. So in principle such things could fill up databases of spam companies who are using robots to search the web for email addresses.

And I wanted to play with cgi a little bit.

"Guests"

2009-10-09T14:04:00.004+01:00

Here is how Argonne National Laboratory, Argonne, IL (just to make sure) treats their "guests":

"The Guest hereby agrees on behalf of the Guest, his/her heirs, administrators,
executors, and assigns to release the Department, Argonne, the UChicago and its Board
of Trustees, and the officers, agents, and employees of the foregoing from all liability for
any personal injuries to the Guest which may arise out of the Guest’s access to Argonne
and the use of its facilities. The Guest further hereby covenants and agrees that neither
the Guest nor the Guest’s heirs, administrators, executors, or assigns will bring any suit
in a court of law or equity or before any administrative agency for any such injury
(including injury causing death) against any of the aforementioned parties."

Did I misunderstand something or is this just crazy and way out of line? The way this is formulated Argonne staff could actively attempt to murder me and I could not sue them. I assume that such a waiver is not in accordance with US law anyway but why do people write such a crap then?

Things like that just piss me off. But I guess it's probably not enough to make me really cancel the trip.

WolframAlpha

2009-08-20T19:00:00.002+01:00

I just got an email from Wolfram|Alpha, particularly about their chemistry interface which is introduced here an here. Wolfram|Alpha is something like an encyclopedia with added functions from Mathematica.

The program can for example easily solve the limit from last post. It seems like a nice free alternative for people who do not want to get a Mathematica license.

It also seems to have enough artificial intelligence to recognize what you want even without strict syntax, for example for this integral.

I am impressed by the mathematics things because I don't know of anywhere else where you can do such symbolic math for free. The chemistry part is also nice because it is all integrated and web based but most of the functionality is also provided by other free programs.

Here is what happens if you enter a sum formula. It nicely presents you the molecular weight and composition.

If you enter a substance like camphor that is in their database than you get a fact sheet. Of course the thing to compare it with, is the wikipedia entry. Most of the numbers are similar. For some reason the NFPA labels are different. I guess the advantage of Wolfram is that it has a more uniform appearance, wikipedia is probably more extensive. And I don't know which one is more accurate.

The tool also does some nice stoichiometry like telling me how many moles are in 5 grams of 95% camphor. There is also a reaction editor, e.g. it balances sodium + water -> hydrogen + NaOH. It even computes a reaction enthalpy but it tells me that the reaction is strongly endothermic. Are just all the signs switched? I mean regular inorganic compounds should not have positive formation enthalpies?

To conclude: Wolfram|Alpha seems to be a nice tool because it combines encyclopedia functionalities with mathematics. The chemistry part also offers good functionalities but nothing really new. If you have any suggestions I am sure they are happy to hear from you at their blog.

Marcus-Levich-Hush equation

2009-07-23T15:07:00.008+01:00

Here is an interesting piece of math:

Or in other words show that in the Marcus-Levich-Hush theory, electron transfer stops when you approach zero temperature. Mathematica tells me that this is the case but I could not really show it. The problem has the following type:

It is not an apparent l'Hospital application and I do not really know what else to do. Anyway, I believe Mathematica and it is also the thing you'd expect physically.

In the full quantum picture you actually have a non-zero transfer rate even when approaching zero temperature that comes from nuclear tunneling (as was derived here and is reviewed here).

Interestingely it is quite easy to find derivations of this quantum formula. But I did not find much about the semi-classical formula, that is shown above. Only that it is some kind of application of Fermi's golden rule which apparently gives this prefactor for the rate equation.

Excited state intramolecular proton transfer

2009-07-18T16:00:00.001+01:00

Actually there is a second article with my name on it that just came out. So I'll use the chance to advertise for it a little bit as well. It is a review (at the time of submission also partially a preview) of excited state intramolecular proton transfer cases. The context is that it is part of a special birthday issue. In fact it is kind of interesting how professors have symposia instead of birthday parties and give each other journal issues for presents ... but that is a different question.

This figure, which is similar to the graphical abstract, kind of shows what's going on: You have a molecule with an intramolecular hydrogen bond (the red atom is typically oxygen, the green one oxygen or nitrogen). You excite with a high frequency UV-photon. Then in an ultrafast transfer process with no or almost no barrier the proton moves to the other side. The new tautomer may either emit a photon with strongly Stokes shifted fluorescence (as for example bipyridyl-diol from the last post) or it may exhibit internal conversion (as e.g. 2-(2'-hydroxyphenyl)-benzothiazole in gas phase).

Bipyridyldiol

2009-07-12T13:15:00.000+01:00

Nice, I just noticed that my Master's thesis article is already available online [1]. Well, you can take a look at it if you are interested in excited state proton transfer or if you want to know what I did for my Master's thesis. I'll put some additional things here that did not make it there. This is basically what happens:

Bipyridyl-diol has two intramolecular hydrogen bonds. You excite it with UV light wait a few femtoseconds and the protons get transferred. It was understood that the double-transfer product DK is finally formed. The main question was wether there were sequential and/or concerted transfers. The general idea was that there would be a branched reaction path: An ultrafast (100 fs) first step that was either a single or double proton transfer and a "very fast" (10 ps) step from MK to DK. According to us it looks more like there is no branched reaction but rather a dynamical equilibrium between MK and DK that cools toward DK. Well I hope some experimental groups are still interested enough in this system to test for this hypothesis.

This is one trajectory, a simulation of the molecule for 300 fs after UV excitation.[2] You can see a very quick initial transfer and then some more transfers.

Actually I wanted also to show the development of the normal modes in the video. To compare them with the results of Stock et al.'s experiment. But this does not seem to work out here because the videos need to have a fixed 4:3 format. So I'll just show a figure. The important thing is that there is strong participation of the totally symmetric modes (blue, red) even if the process does not conserve the symmetry. Another very interesting thing is that activation of the non-totally symmetric (black) mode is a violation of the Franck-Condon rules. A way to explain this is that the Franck-Condon rules work only under ideal assumptions and not with a strongly anharmonic reactive potential.

Here is another trajectory for comparison. In this case the second proton transfer occured only a little bit later.

Actually another nice figure would be this one. What I am doing is projecting the trajectories onto a normal mode. And then I can average for every time step over the 36 trajectories that we ran. This time-dependent average should represent the coherent motions. Here I am showing 17a_g, an aromatic breathing vibration, which is the classical case for a coherent Franck-Condon excitation (in the context of proton transfer the lower frequency skeletal modes were of more interest). In the harmonic vibrational analysis that we did at the DK equilibrium geometry, the mode has a wavenumber of 682/cm. This corresponds to a period of about 49 fs. Well and there really is a coherent oscillation with just that frequency. So we see that the harmonic vibrational analysis at the minimum and the dynamics nicely work together. If I compute the standard deviation over time of this time-dependent average then I get only one number per normal mode. These numbers are what we are showing in Fig. 10. And by the way: The tools to do this are in the new Newton-X version (aside from many other nice things ...).

[1] And interestingly there is a direct link to facebook which of course I had to click.

[2] One of these 300 fs RI-CC2/SVP-SV trajectories takes about a month on one processor.

Excited state H transfer

2009-07-03T16:39:00.000+01:00

I liked the introduction of this article by Röthlisberger because it nicely explains the processes in excited state proton or hydrogen atom transfer. Most of it is well explained in this Figure.[1]

In the ground state the n and π orbitals (shown on the left-top and right-bottom) are each doubly occupied. Excitations into the two virtual orbitals shown (left-bottom and right-top) lead to three states of interest: ππ*, nπ*, πσ*.

In the ππ* and nπ* states it can be seen that electron density is shifted from the O to the N. This increases both the acidity of the O and the basicity of the N. In the cluster shown this induces proton transfer through the ammonia molecules. Actually there is a very nice movie showing this transfer in their supporting information.

The situation is completely different in the πσ* state. If the anti-bonding σ* orbital is populated, the bond is no longer stable. The molecule stabilizes by dissociation of a hydrogen radical (i.e. hydrogen atom). If the hydrogen atom takes part in a hydrogen bond, you can have excited state hydrogen atom transfer. The orbital corresponding to this is shown in the left bottom. It is very diffuse and has probably also some Rydberg character.

When I decided to write the post, I had only read the introduction which is very nice and helped me finally understand the difference between excited state proton and hydrogen atom transfer. The part that seems kind of strange is that they only computed one trajectory. And for this one trajectory they had 1024 processors on a Blue Gene/L. With atom centered basis sets in Turbomole you could almost do it in real time if you had 1024 CPUs (or pretty fast anyway).[2] So it may be a problem with the plane waves. Or I got something wrong. Or I am just jealous because I never had 1024 CPUs at my service.

[1] Call it advertisement ...

[2] Not quite real time as one CPU cycle is already about 3 orders of magnitude longer than the process observed.

Randomness

2009-06-26T21:00:00.001+01:00

If you wonder why you lost all your money in the financial crisis or why people get it wrong when it comes to statistics, here is what could read: Fooled by Randomness and The Black Swan by Nassim Nicholas Taleb. The first one is about how we are unable to understand chance intuitively and the second one particularly about the rare event that may compensate everything that came before it. And he has a really nice writing style.

According to the "narrative fallacy" you are more likely to believe my point if I tell you a story, alright: Let's say you are playing "Mensch ärgere dich nicht" or some other game with a die. For many rounds now you have not been able to get a 6. Are you more likely to get a 6 in the next round? "Next time it really has to be a 6!" But of course the chance is still 1/6. Or to apply Taleb's street smarts logic: Someone is probably cheating and the chance for a 6 is even lower.

To stay with games, there is also what Taleb calls the "ludic fallacy": People tend to think that the kind of randomness that matters is the one that follows a known statistical distribution. But it is the events that no-one thought about that matter, or the ones that do not happen in "ten billion times the age of the universe according to our model" ... the Black Swans.

So what is the problem: of course evolution. Our rational thinking (or "Reflective System") is fairly new in evolutionary terms and it does not really play a role when it comes to making decisions (at least not as much we would think it does). What really drives us is intuition (or the "Automatic System"). And statistics in the Automatic System come in the sense of heuristics.[1] The problem is that these heuristics only really apply to Stone Age conditions. For example they don't take into account that the person on whose clothes you just spilled your coffee cup is a complete stranger that will be 5000 miles across the ocean within a few hours. So you will apologize and try to help even though this behavior is probably nicer than what would be purely rational. In this case it was good but the problem is that our intuitive heuristics fail in many cases.

Actually this brings me to another nice point which is interesting, at least from a slightly academic point of view. A major question in trading is if it is possible to be on average better than someone who buys purely at random.[2] Let's assume some stock is lower in value than it should be according to some economic considerations. If there is a way to find out, people will buy the stock and the price will go up until the stock is not cheap anymore. Therefore it is very difficult to get rich in the stock market through information unless you are either very smart or very quick. But in our real world stock prices are not ideal but biased by humanness (e.g. because of people that think that they are very smart or very quick, or generally herding phenomena, and so on). So stock prices will reflect the real value with a bias by human (Stone Age) intuition. All you have to do to get rich, is counter-intuitive trades. I like the thought but the problem in economics is always that oversimplified theories lead to elegant results but to the destruction of trillion dollar hedge funds. Anyway in essence my point is similar to what Taleb says only that he gives the examples.

[1] Some of these things are nicely explained in the first part of "Nudge" by Thaler and Sunstein. In the second part they will probably explain how to nudge people to do the right thing but I am not there yet. "Nudge" is also a nice book but it is a bit too academic to be an easy read ...

[2] Here the point is average. You can (and should) affect the deviation by diversification. In case we are talking about your funds you should keep the deviation low. If it is someone else's money, you should gamble. If it goes up use the hype to get rich. And when it crashes, you should disappear or take the settlement, depending on your contract.

Antisymmetry

2009-05-31T17:00:00.002+01:00

Antisymmetry, second quantization and things related to it are something I kind of stayed away so far. And now that I have finally taken a look at it, I notice how nice the math behind it actually is. And the important things that many things work in general and do not require orbitals.

We consider functions of the following form where A is some set and N a natural number[1]:

Let V be the vector space of all such functions. And let's look at operators [2] acting on V. First we could consider the Transposition operator

This is a well defined mapping and the application is straight forward, e.g.

And if you think about it some more you notice that it is linear.

And if I am correct it is both Hermitian and orthogonal with the usual skalar product. This would mean that all the eigenvalues are either 1 or -1. Either way, we are interested in the eigenspaces corresponding to the eigenvalue -1, i.e. functions that are antisymmetric with respect to the transposition. The intersection of all these eigenspaces V_A is the set of functions that are antisymmetric with respect to all the transpositions. It is a vector space since it is formed as a intersection of vector spaces.

All fermion wave functions that comply with the Pauli principle have to be taken out of this space.

Next it helps to be a little bit more general and to introduce the permuation operator (related to a permuation σ out of the symmetric group S_N) as a generalization of the transposition operator

This operator is also linear and I think also Hermitian and orthogonal.

Then you could write the antisymmetric space also like this

since every transposition is a permutation with negative sign and every permutation can be formed out of transpositions and the sign corresponds to how many you take.

Well let's define one more operator, the antisymmetrizer

You can show that this operator is a projection operator. This is equivalent to the condition that applying it twice is the same as applying it once, since the vector is already projected after the first time. The proof is also straight forward. Applying it twice leads to a double sum

But this just means that you sum n! times over all permutations, and it turns out

And next you could prove that the space it projects into, is just the V_A we had before. You have to show that any projected function fullfils the condition for V_A above and that every function which fullfills the condition remains unchanged which is both quite straight forward.

Maybe I should also add why these results are nice. I did have the expression for the Slater determinant before. But by considering that this is nothing but a (normalized) projection of an orbital product into the antisymmetric space it is easier to deal with it. It is no longer some weird expression but nothing but a linear operator. And things like adding Slater determinants are much clearer when you consider this. And that is why it is cool.

And actually as a next step you can consider commutation. These operators commute with the Hamiltonian because they are nothing but relabelling of equivalently treated electrons. Only for that reason you know that eingenfunctions of the Hamiltonian which are also antisymmetric even exist. You restrict the Hamiltonian eigenfunctions to the ones with a -1 eigenvalue of the antisymmetrizer. In the next step you may also do some restrictions according to spatial symmetry - the interacting spaces will transform like irreducible representations of the symmetry group. And finally spin. In the spin case the interacting spaces transform like representations of the unitary group (and I am kind of trying to understand why that is).

[1] Physically spoken N is the number of particles and A is the set of possible coordinates of the particle.

[2] As I probably said before: An operator is a function that makes a function out of a function. Or if you don't like the word "function", it is a mapping between two vector spaces.

Hubble

2009-05-14T09:47:00.010+01:00

For some reason CNN is one of the only channels that I can access from my room here in Prague that is not either in czech or featuring arabic telephone sex commercials. One thing I saw there were some pretty amazing bulls and bears that reminded me of N. N. Taleb. But I also saw some exciting space shuttle features. I like space shuttles but I am wondering what they are actually good for. (I am not an expert on space shuttles but I did once hear an astronomer talk about them.)

As I understand it, a space shuttle is like a Porsche. It's exciting, it's prestigious, but it does not really get you anywhere you could not go without it. (The advantage is that a space shuttle is not as noisy and does not cause traffic jams that I have to carefully bypass on my bicycle without collecting any mirrors.) The question that the people at CNN never asked any of the space shuttle experts is how sending up a 7 austronauts and assembling a telescope in deep space compares to assembling the telescope down here and sending it up by itself. I guess a space shuttle is a nice piece of science fiction, without the need of fiction, but probably not the most cost efficient way of having a telescope in the sky.

CERN

2009-05-12T12:25:00.006+01:00

I should not critize my home country but this was a masterpiece: We are apparently quitting CERN! Especially now that they build the LHC.

Well the upside is: When the black hole is created that swallows the whole world, we can say it was not our fault.

Edit: and the second upside is that they will have about 20 million extra euros that they can put into my stipend.

Edit2: We are not leaving CERN thanks to public pressure.

DNA

2009-05-11T14:21:00.005+01:00

Ever since Jurassic Park we know that DNA looks pretty cool. Here is a piece of B-DNA.[1]

Or with thicker sticks:

In fact there is of course not so much space between the bases. Here is what it looks like when you draw the van der Waals spheres for the bases. It is actually quite difficult to model this stacking interaction because it is almost pure dispersion. From a computational point of view, dispersion is electron correlation. If the effect you are interested in is pure electron correlation, then you'd better model it really well. Hartree-Fock gives you zero dispersion. MP2 is great if you kind of want to include some electron correlation but it is not accurate enough in this case. What they are actually doing is high level coupled cluster CCSD(T) extrapolated to the complete basis set limit.

Another interesting question is how defects are propagated in such a framework: holes, electrons, or electron hole pairs. First of course you want to preserve the integrity of your genome. And second you want to make nano-robots. DNA is already a self-assembling structure with molecular recognition. If charge transport is better understood and conductive analogues are found, DNA will kick nanotubes' ass.

By the way: for printing out pymol graphics (see again [1]), the ray_trace_mode setting is nice. Especially:

set ray_trace_mode=1

Then you get some nice black frames instead of fuzzi ends in the print.

[1] from the crystal structure 1BNA. And drawn with pymol. Maybe I should eventually switch from pymol because they are becoming more and more commercial. But then I guess they would not mind me as a little personal user.

Check every day

2009-05-07T17:21:00.005+01:00

I have three more items that you should put on your "Check every day before you start working list":

What should already be on there:

exp(K)

2009-05-05T16:52:00.008+01:00

It is time for some nice extra math, with some proofs and applications. I have two theorems that don't look obvious at first. And initially I kind of tried to avoid them.

The first one is that the matrix exponential of a real antisymmetric matrix is an orthogonal matrix with determinant +1

The second one is that for any matrix K (for example the anti-symmetric one from I) and a matrix A the following holds.

Or explicitely:

This is an expansion in the powers of K. And if K is small, you can cut it at some point. The reason why this is useful is that it gives you a direct possibility to manipulate orthogonal matrices. You don't have to worry about making the matrix orthogonal but you just choose any paramters for K and exp(K) will be orthogonal. II gives the possibility to carry out the basis transformation in an efficient way.

The first part of I is actually very straight forward if you write it down like this. It is provided that the limit of a transposed sequence is the transposition of the original limit and that the inverse is formed as shown.

I will prove the second part with the spectral theorem (maybe you can also do this in a more direct way). If we move into the complex space. K is now correctly spoken anti-Hermitian, and we know that a unitary matrix U exists which diagonalizes K. The eigenvalues of K are purely imaginary (yet quite important).

The exponential is given according to

At this point you can also notice that exp(K) is unitary because it is formed as a product of three unitary matrices.

Can we say anything about the determinant? Yes of course.

The determinants of the mutually inverse matrices UH and U cancel out. The determinant of a diagonal matrix is just the product of the diagonal terms.

The sum of the eigenvalues is the trace of K. In a real antisymmetric matrix all diagonal elements have to be zero (therefore also the trace).

And that is what we wanted to see.

The proof shows that exp(.) is a mapping between the set of real antisymmetric matrices and the set of real orthogonal matrices with determinant +1 (and of course that -exp(.) is the mapping to matrices with det=-1). In fact it would also be necessary to show that this mapping is bijective (i.e. that there is a 1:1 correspondence). This is apparently true but I don't think the proof is so simple. Interestingly it does not work with complex anti-Hermitian matrices. In this case there are n extra purely imaginary values in the diagonal. And the mapping is no longer injective for example:

Point I was the elegant part. Point II is more of the index juggling variety but a nice piece of complete induction. I will probably put that up soon, too.

The Singularity is Near

2009-04-29T15:30:00.001+01:00

If you want to know what's going on in the future you could read the book "The Singularity is Near" by Ray Kurzweil or check out his homepage. The Singularity is the point in time when artificial intelligence passes by human intelligence in all respects.[1] According to Kurzweil this happens about 2050. By then it is literally unnecessary to think. The process is probably gradual. We have long expanded our memory capacities by writing things down. Computational capacity is of course expanded and there are many things which we just don't learn any more, like taking square roots. Computers are already better chess players. Eventually computers will be better at everything but we will probably stay close connected.

Kurzweil describes three revolutions: genetics, nanotechnology, and robotics. You can first extend your health with genetics and biotechnology. The author takes something like 120 dietary supplements and 7 IV's a week. And because of that he did not really grow older during the last 20 years. Later you can use nanobots to fix everything that may have gone wrong. Finally you can upload your brain to a computer. And if you take frequent backups of your data, you will basically remain forever.

If you actually uploaded yourself, then of course your virtual appearance will be everything there is. Maybe that will be called web 3.0 or web 4.0 . But it cannot be denied that the web is moving into that direction.

The idea of wikipedia was to give everyone free access to the sum of all human knowledge. And actually it became all human knowledge in a matter of seconds to everyone at any place. Thanks to search engines using artificial intelligence and mobile internet, respectively. Just an example: I was arguing with my brother about the sugar content of tonic water. Ten years ago we would have gone to the store, bought a bottle, and looked it up. Now we could hear the shocking truth (tonic water has the same amount of sugar as coke) in a matter of seconds while driving on the high way. What will it be like ten years from now? The change will probably be even more drastic thanks to the law of accelerating returns. Maybe something like a mobile broadband brain tap internet connection.

So how will we react to artificial intelligence? I guess we will respect its feelings. Actually it even feels weird to say something mean to the IKEA lady. So we will probably be nice to a well programmed AI. The question if it is conscious is something different. The problem is that consciousness is something that cannot be scientifically adressed. Or how would you prove to someone that you are conscious?

So what will happen after that singularity? We will send nanobots around space. They will either go close to the speed of light. Or if the computers are smart enough, faster - e.g. through wormholes. And since there are no nanobots buzzing around everywhere, there are probably no aliens around (unless they chose to remain hidden).

Anyway, the point is that the future is probably totally different than it is now. And it probably makes sense to think more about the "deeply intertwined perils and promises" that these developments will bring. Technology will probably keep improving people's lifes. But weapons will also become more powerful and more difficult to control in case something goes wrong.

[1] Maybe you could also look at accelerating paradigm shifts. If they form a convergent series then there is one point in time where an infinite number of paradigm shifts will have happened. (Or at least quite a few if the convergent model is not quite correct.)

No Fluorescence (2)

2009-04-08T20:15:00.000+01:00

So how is non-radiative decay possible? I showed the conical intersections of cytosine some time ago. These geometries, where the first excited and ground state are of the same energy, exist. But the chance of exactly reaching one is zero (since they are only an N-2 dimensional hyperline in the N-dimensional space of geometries of a molecule).

The point is that the picture of isolated electronic states breaks down, i.e. the Born-Oppenheimer approximation. In solid state chemistry that is what they call electron-phonon-coupling.

If you consider the nuclei in terms of an external perturbation, you can take at look at the adiabatic theorem: If nuclear motion is fast and there are close lying electronic states, the states will mix. You get the fast nuclear motion from the excitation energy. About close lying states: Typically all the excited states are close and the molecule quickly reaches S₁ if it was initially excited to a higher state (Kasha's rule). If there are suitable intersections to the ground state it can even relax completely.

For doing molecular dynamics you need the energy of your electronic state Ψ_i as always. But you also need couplings to the other states:

They are called the non-adiabatic coupling vectors and give the major post Born-Oppenheimer contributions.

Here is a little dynamics movie. This is a test MCSCF run on cytosine. The electronic structure is from the Columbus program and the non-adiabatic surface-hopping dynamics are with Newton-X (both from our group actually).

The idea of surface-hopping is that we still want to have a classical trajectory as a basis, but post Born-Oppenheimer corrections are introduced through jumping in between states. Here you can see the cytosine molecule with color coded electronic states: green S₁, orange S₂. You see how the molecule jumps into the second excited state intermediately and then stays in the first excited state. There is a lot of motion because of the high excess energy.

There is no decay to the ground state in these 341 fs simulated. It should happen soon after or tanning would be extremely dangerous. But the time simulated is still rather short, it is for example a ten-thousandth of the typical fluorescence life-time.

Symmetry elements

2009-03-22T13:15:00.007+01:00

Here is another awesome tool that I just saw in a Chemical Forums post. It lets you visualize symmetry elements and even animate the corresponding symmetry operations.

Let's take for example trisoxalato-iron(III). What symmetry elements does it have?

First there is a C₃ axis that goes between the ligands. It goes through a face of the coordination octahedron.

Second there is a C₂' axis perpendicular to the main axis. It goes through an edge of the octahedron.

There are no more symmetry elements, the complex is of D₃ symmetry.

Druglikeness

2009-03-08T13:08:00.004+01:00

Because I like cheminformatics tools and because I got a nice e-email, here is one: Osiris Property Explorer. It lets you predict pharmaceutically relevent properties of a molecule and it computes absolute configurations in chiral molecules.

This is what it looks like with my test molecule.

You can see that the two chiral centers are of S configuration.

More important it has a druglikeness of 1.8 and no Toxicity Risks. So it would make sense to pay me huge amounts of money for a patent and start building a plant for large scale production. Because with this kind of a druglikeness it does not matter much that it has a high octanol/water partition coefficient, is not soluble in water, and has a rather high mass ...

Divalent hydrogen

2009-02-15T17:28:00.009+01:00

In a recent Nature communication scientists from Harvard presented the preparation of divalent hydrogen. It goes according to the following scheme:

Tell me I am picky and jealous. But what should you expect of this? Is only this figure the only thing wrong in the article? Take it as a mild critic but I am disappointed. There is probably great science behind all that and mistakes happen but should you not proof read before you publish in a high impact journal? And it really takes away information because it is not obvious if the imine or the amine is present.

Maybe the correct way would be to silently tell the author. But if it is a Nature paper it should withstand some public critic in my eyes.

Since I am already complaining: Not stopping at a pedestrian crossing is like aiming a loaded gun at un-armed people who are already abiding to rules that only have to exist because of the car drivers.

Fluorescence (no)

2009-01-17T21:30:00.004+01:00

For the longest time I did not understand how it would be possible that a molecule should not show fluorescence. You take a photon to push it up to a higher level. So how should it be able to move down in a different way, between two different quantized levels. The main errors in this thought are the static view point and the Born-Oppenheimer approximation.

The first thing to consider is that gaps between electronic states (in the Born-Oppenheimer approximation) are not the same at every geometry. At the ground state minimum there is typically a gap of several eV to the first excited state and correspondingly excitation is in the UV. If you distort the geometry, the ground state energy increases (since we were at the minimum). Since the minimum of the excited state does in general not coincide with the minimum of the ground state, it can relax to a lower excited state energy. Ground state and excited state move together and the gap decreases. Fluorescence will be at a lower energy than absorption. This is just the semi-classical explanation of the Stokes shift.

You can go on trying to decrease the gap. And it turns out that for many molecules of interest there are actually geometries where the ground and first excited states are degenerate. These are called conical intersections. And there is not just one such geometry but in the general case the crossing seam forms a hyperline in the space of all nuclear geometries (i.e. an N-2 dimensional subspace, where N is the number of nuclear degrees of freedom). Conical intersections are usually structures strongly distorted from the ground state geometry (hence the ground state energy increases). And they are adapted to the excited state (to keep its energy low).

Nucleic acid bases are of interest in this context. As I read it in a nice popular science article: When a photon from the sun is absorbed by a molecule, this essentially means that this molecule is heated to the temperature of the surface of the sun. Conjugated polar molecules like nucleic acid bases definitely absorb UV light. So why are our cell nuclei not scorching away when we are out tanning a little bit? Well the answer is that they are very efficient in giving away the excess energy and thereby returning to the ground state.

Actually people are arguing that photostability was one main points of selection pressure in the early biosphere. A support of this is that the nucleic acid bases have very short decay times compared to analogues with different substitution patterns. Analogues were destroyed by photochemistry while the bases that are now actually in use could protect themselves with photophysics (i.e. non-radiative decay).

Cytosine is an interesting example [1]. In this case there are three excited states of interest and all of them cross with the ground state. The question is which one is actually the decay channel and there is not really consensus. Well I just want to show the geometries.

The ground state minimum is almost planar with same pyramidalization on the amino group.

This is one of the conical intersections, usally called "twist CI". The excited state is a biradical state localised mostly on the two C atoms in front (coming from the ππ* state). Through this twisting the biradical state is stabilized because the electrons move out of each other's way. I think you could also say that you have pyramidalized radicals just as usual in organic chemistry.

This on is the "sofa CI" which is related to the n_Nπ* state (the n orbital localized on the N that is not protonated). Ring puckering and stretching of the amino group stabilizes this state.

The last minimum of the crossing seam between first excited and ground state is this one. It comes from the n_Oπ* state.

Well these conical intersections exist. But the probability of actually reaching one is zero. As I said crossing seams are hyperlines. Imagine a line in the room where you are just sitting - any given random point will at least have a little bit of distance [2] from the line. The way out is that with close lying electronic states the Born-Oppenheimer approximation breaks down and a transfer between electronic states is possible through coupling with the nuclei.

[1] The method was MCSCF in Columbus but I cannot really give computational details here. This is more thought as a "community outreach" - just showing some nice pictures.
And you can for example check out this ref if you want to know more about cytosine.

[2] an ε > 0, if you like calculus

Proton Transfer

2009-01-15T21:40:00.006+01:00

This time I will show a few "quantum images". Hydroxybenzoquinoline is kind of a cool molecule because it shows excited state intramolecular proton transfer. After UV excitation the proton is transferred from the oxygen to the nitrogen. This happens on a time scale of about 30 fs because the process is essentially barrierless. This is a time scale where you can do ab-initio molecular dynamics. And with modern LASER pulses you can also examine it experimentally. That's why these molecules have been rather popular recently.

What I want to show here are the orbitals because I think they look kind of cool. What you see are the two highest occupied orbitals (n, π) and the lowest unoccupied orbital (π*). The two important excitations are just between these orbitals leading to the ππ* and nπ* states respectively.

Related to the proton transfer it can be seen that the π* LUMO has more density on the N which increases its basicity to help it catch the proton.

Something else that is interesting is that a non-bonding orbital (coming from MO theory) really looks a lot like the free electron pair you would talk about in VB theory. I was arguing with an analytical chemistry professor once that he should not interchange the terms non-bonding orbital and free electron pair because they come from different approaches. Now I see that he was right but I did not believe him then because he was from analytical chemistry.

Enol (ground state)		Keto (ππ* state)
	π*
	π
	n

Just a short summary of the software: the orbitals are actually from semi-empirical DFTB. And the pictures are made with VMD.

Normal modes (3)

2009-01-01T15:40:00.003+01:00

Ok finally comes the part that I wanted to show. But first a little jmol example to provide you with something to play with.

What you see are the normal modes and frequencies determined as I will describe below. They are pretty much what you would expect from IR spectroscopy. Only that they are a little bit higher which is mostly due to the fact that the harmonic approximation is not quite correct. [1]

You can rotate the molecule by dragging with the left mouse button and switch between the vibrations with the right mouse button and changing the "model".

What we want to do is solve this coupled system of differential equations:

The idea is that if we bring it do diagonal form, the equations are uncoupled. So what can we do? One way would be

But the problem now is that M^-1H is unlike H not a Hermitian matrix and we do not know if it can be diagonalized. But you can rewrite it as:

Or if you introduce (square root) mass weighted coordinates S(t) it looks like this.

H_S is actually the Hessian matrix in mass-weighted coordinates. The important point is that it is Hermitian and hence an orthogonal basis of its eigenvectors exists, i.e. there is an orthogonal matrix Q and diagonal matrix D with the property

Plug this in, multiply with Q and that's actually the result.

The idea is now that you have uncoupled collective coordinates and you can solve each line of this equation system separately and get some cos function each like in the first post. The eigenvectors of H_S (the columns of Q) are the corresponding motions. And the eigenvalues determine the frequencies.

[1] The lowest vibration is kind of out of place. With DFT/B3LYP cytosine is planar at its minimum geometry and the NH₂ pyramidalization has a very weak force constant. Actually with MP2 and some other ab-initio methods the amino group is pyramidalized at the minimum. So I would guess that DFT is wrong. It is not a big issue but you can see that DFT is sometimes a problem and it makes sense to check the results a little bit.

Normal modes (2)

2008-12-13T13:22:00.005+01:00

Now I am used to the fact that almost everything in theoretical chemistry comes down to eigenvalues and -vectors. But it kind of struck me when I first heard that about normal modes some time ago. In short: normal mode frequencies are the eigenvalues and normal modes the eigenvectors of the Hessian matrix of the energy in mass-weighted coordinates.

Here's the math behind it because it's always fun to make some formulas in LaTex.

We look at the effective potential energy V(R) for the nuclei of a molecule with N atoms in the Born-Oppenheimer picture [1] which is a function of the nuclear positions R=(x₁, y₂, ..., z_N)^T. This function is expanded into a Taylor-Series up to second order which you can write like this (where a suitable origin 0 is chosen):

The first term is called the gradient, defined according to

The second term, the Hessian matrix

That was actually just kind of a warm up as we are interested in the energy gradient which we can expand to first order in the following way:

With normal mode analysis we can describe the motions at a local minimum. At a local minimum the gradient is 0 and all the eigenvalues of the Hessian are greater or equal to zero. The first condition gives:

Now we have a convenient description of the gradient that we can plug into Newton's second axiom. We get a differential equation system with an equation for every coordinate. This can of course be represented by a matrix equation. First it is convenient to introduce a diagonal matrix M that contains the masses (each of them 3 times because there are x,y, and z coordinates for every atom):

Then the equation system looks like this:

If you write out the matrix equations then you have a system of coupled ordinary differential equations. And actually the part that I wanted to show is how you can uncouple this system by taking proper care of the mass and diagonalizing the Hessian. Well next time ...

[1] That means electronic kinetic energy and all potential energy terms.

Normal modes

2008-12-08T12:00:00.002+01:00

I am kind of excited because I finally understand how normal mode analysis works. Finding the normal modes and corresponding frequencies (and intensities) of a molecule is what you need to predict or interpret IR/Raman data. Like so many other things it comes down to eigenvalues and eigenvectors.

We start out with Newton's second law

The first consideration is the one dimensional oscillator. We have one variable x, and the force defined according to:

This leads us to

In other words we are looking for a function whose second derivative is the same function with a minus ...
How about sine or cosine?

The following satisfies this equation. The proof is just differentiating it twice.

It is not quite so obvious that with two real numbers A and δ this is also the complete real solution of the differential equation (where A is the amplitude and δ the phase) but also not so important in this case.

The part that has an immediate application is which gives the frequency. k is the force constant and it increases according to: torsion < bend < stretch (single < double < triple). m is the mass. The highest frequencies are for low mass and high force constant, i.e. X-H stretch. Triple bond stretches are high because of high k, then double bond stretches and so on.

This was the one-dimensional case. Next time I want to show how to reduce the general case to isolated one-dimensional equations like shown here. And that's where the Linear Algebra comes in.

Ribonuclease A

2008-11-24T17:00:00.002+01:00

Here is some stuff from my last lab [1]. The objective was to do force field calculations (with CHARMM) on our protein of choice. What I picked is Ribonuclease A, the enzyme that digests RNA in our food.

Ribonoclease A is a sturdy little molecule with 124 amino acids. Its four disulfide bridges (yellow in the picture) make it very stable and comparably easy to isolate. Just 124 amino acids also makes it nice for computations. The pictures were made with VMD which makes nice pictures and is also a pretty versatile tool that reads tons of file formats.

The active group of Ribonuclease A comes from two Histidine (black) residues that help hydrolizing the phosphate in the nucleic acid. Lysine (orange) provides the positive charge to interact with negatively charged phosphate groups.

The first thing you can do with a molecular structure is an energy minimization. These are the C_α traces of the crystal structure (yellow) and the vacuum optimized structure (red). They are in fact quite different even though the red structure is just the next local minimum.

The major problem with this optimization is that I did not include solvent effects. The electrostatic interaction energy (in the macroscopic limit) is given by

The crucial part is ε_r, the relative dielectric constant of the medium. In water ε_r=78. In other words: charges are strongly shielded compared to a vacuum. In the microscopic world things aren't as simple but it is still the same trend. The difference is probably a major reason why the structures are so different.

Another thing you can do with structures is normal mode analysis. That means computing the eigenvalues/eigenvectors of the mass weighted energy Hessian matrix which correspond to the frequencies/motions of the normal modes. Ribonuclease A is small for a protein but it still has 1856 atoms. So we have 3x1856=5568 degrees of freedom. We have to compute, store and process about 30 million matrix elements. This seems much but these days most regular computers will be able to do this in less than half an hour. 30 million matrix elements corresponds to something like 300 MByte, nothing too crazy. [2]

These are the three lowest frequency modes (aside from translation and rotation) of Ribonuclease A. They are large scale backbone motions. A weak force constant and a high effective mass causes them to have wave numbers below 10/cm.

[1] Not the kind of lab where you need a lab coat.

[2] Actually if you are just interested in a few eigenvalues and -vectors, you can go much higher. For example multi-reference configuration interaction with Columbus (which is a great program by the way....): finding the first 4 states out of a 8569890 x 8569890 matrix was possible within about 17 hours. In this case you don't store trillions (7E13) of matrix elements and do a diagonalization. You use the iterative Davidson procedure which will find an eigenvector close to a starting guess.

GFP pictures

2008-11-09T22:30:00.000+01:00

Time for some more graphics. I picked GFP which is not really a new idea but it is a pretty cool molecule.

In the sphere model all proteins pretty much look alike. But I still think kind of cool, especially if you use Lightnir's QuteMol preset for pymol.
The cartoon model shows you more about the structure. GFP consists of a β-sheet "barrel" that encapsulates the chromophore. You need a rigid structure if you want to make sure that the chromophore does not quench before it is able to show fluorescence. It is interesting to compare this to Rhodopsin which is apparently much more flexible. This makes sense since in Rhodopsin the biological activity is related to non-radiative decay and isomerisation.

The chromophore is trapped at the center of the barrel.

Interestingly the chromophore is made directly out of amino acids. You can take a look at it and think of how this is done.

The answer is here. The hydroxyphenyl comes from tyrosine (as expected). The imidazolon ring is formed after cyclization, the bridging double bond through oxidation. It is interesting to consider that this apparently works in many different organisms and not just in the jelly fish aequoria victoria. The driving force seems to be that the amino acids are pressed together in the barrel.

The chromophore is part of the protein chain. It is connected to the helix that goes through the center of the protein.

Or from a different perspective.

QMC@home

2008-11-04T21:03:00.003+01:00

What you need in computational chemistry is - not so very surprisingly - computers. You can either use your own little cluster. You can try to get time at a super computer. Or you can have other people do the work for you.

Actually I am advertising for the competition a little bit but I think the idea is pretty cool. You are probably used to distributed computing with the SETI project. QMC@home does the same in quantum chemistry. The disadvantage of their Quantum Monte Carlo approach is that you need a huge amount of computer time. The advantage is that you can easily distribute it over many computers with little communication between them. The ideal case for voluntary grid computing. It is organised over the same BOINC platform that is also used by the SETI project. So if you download it and get bored by Quantum Monte Carlo, you can switch the system and start looking for extra terrestrial life.

QMC@home has this screensaver that tells you what's going on. (In principle it should have it, it's not working on my computer).

In this case they are computing interaction energies for the cytosine dimer. Van-der-Waals dispersion interactions are a difficult problem because you need a highly correlated wave function to describe them correctly. Maybe QMC is a good way to do it.

I like the idea because it seems like a waste to have so many unused computers standing around everywhere. It has to be said though, that the computation uses up energy. But then when I look at people that don't even turn off their screens when they leave their computers, it does not seem like they would complain about a little bit of extra energy use.
It would probably be a good idea if they had some kind of reimbursement system. But on the other hand it is kind of fun as it is and you get a sense of accomplishment just for keeping your computer on for some time. And now I also have an excuse for not turning it off every time I leave my workplace.

Iodine dynamics

2008-10-27T14:00:00.000+01:00

I summarized the math in the last post. Here is what you can do with it.

I₂ can be described with Morse potentials. Here we are looking at the electronic ground state and the second excited state (because that is the bright state). Energy is plotted against I-I distance. I guess it is pretty much straight forward to solve the stationary Schrödinger equation of such a system to get the vibrational levels and wave functions. The lowest 20 ground state levels and 80 excited state levels are shown.
This is the textbook situation of a larger equilibrium bond distance and lower binding energy in the excited state. Vertical excitation from the ground state will put the wave packet on the left edge of the potential in a highly excited vibrational state. In the classical picture you can say the molecule will start vibrating because the two atoms are closer together than the new equilibrium distance.

Excitation from the v=0 vibronic ground state would require about 509 nm excitation energy and would lead to an almost dissociative state. In the simulation we started with v=1 where you have more density at a larger bond distance and a vertical excitation of 588 nm from the second maximum. v=1 is a realistic situation because this level is just 213/cm or 2.5 kJ/mol above the ground state, so we have 36% in this state (relative to 100% v=0). The situation works together with what we observe. 588 nm corresponds to yellow light. Yellow absorption gives a purple appearance.

In the simulation there was a 300 fs Gauss pulse (simulated with numerical Runge-Kutta integration) and then the wave packet evolved like I explained in the last post. With this ultrashort pulse, excitation of a coherent wave packet is possible. This wave packet is formed above the second maximum of the ground state function, then it oscillates back and forth. These oscillations occur with a period of 333 fs on average, this fits very well with the energy gap between the two levels where the wave packet should be mostly localized according to the excitation energy.

The classical explanation for the oscillations is that since the pulse is only about the length of a vibrational period or shorter, you can excite coherent motion.

You get the same with quantum mechanics but with a different explanation. First you look at the energy uncertainty of the pulse:

Another formula on Wikipedia looks almost the same, the energy of a harmonic oscillator:
In other words: The energy uncertainty of a pulse corresponds to the energy gap between two levels of a system with an oscillatory period equal to the pulse length. If the pulse is much longer than the oscillatory period, excitation will be sharp - one eigenstate will be excited and all subsequent motion is only a phase factor in the complex plane. If it is the same length or shorter there will be excitation into several levels and temporal evolution like I explained last time.

The software used is: Fortran scripts for the numerical integration, the Python Pylab package for creation of graphics, Video Mach for making a movie out of the pictures.

TDSE

2008-10-25T15:29:00.007+01:00

In a way it is easier to solve the time-dependent Schrödinger equation than the time-independent one. If you don't look too closely at H, it is a first order ordinary differential equation. And not even an eigenvalue problem.

If H is constant over time the solution is:

You are probably used to this if you have multiplication with a constant. It works the same way with a linear operator like it is the case here.

This way you can propagate a wave packet when you know the potential energy and the wavefunction at the start.

For every Hermitian operator (like H) there is a basis of its eigenfunctions. The physical interpretation is that the state of the system can always be seen as superposition of eigenstates.
Expand Ψ into these functions. And apply the fact that they are eigenfunctions.

You have time-independent functions with constant weights, only their phases change. The wave packet is the interference pattern.

The harmonic oscillator is an interesting example. The energies are E_k = (k+1/2)hν. If you plug this in, you will notice that it is a periodic function and the frequency is ν the classical frequency of the oscillator and the frequency of the light that causes the transition. (To be precise: Ψ*Ψ is the same after 2π/ν and Ψ has the opposite sign, it has a period of 4π/ν.) I will show in the next post what such a wave packet looks like.

One more consideration. What happens to an eigenfunction in the time-dependent formulation?

It shows a phase change but the physically relevant quantity Ψ*Ψ stays the same. We have a stationary solution.

Are all enantiomers equal?

2008-10-02T20:09:00.004+01:00

The typical statement is that the physical properties of two enantiomers are exactly the same. Differences are only in interactions with other chiral substances and polarized radiation. This is based on the intuitive assumption that a mirror image of our world would behave exactly the same as our world. In physics this is called parity symmetry [1]. But interestingly parity symmetry is not always conserved.

Out of the four forces of nature, it's the weak one that steps out of the line and breaks P- (and even CP-) symmetry. As the name suggests it is weak. And it is also short range, on the order of attometers. That's why we don't notice it. Nonetheless theory predicts energy level splittings between enantiomers. This could be on the order of 10^-15 cm^-1. This is apparently extremely small and even the newest experiments haven't gotten below 10^-13 cm^-1.

Don't bet too much money on the emergence of weak force intermediated enantiospecific synthesis yet. But even though P-violation is out of every day life it seems that with improved experiments or in different systems the effect could be observed. And then chemistry could be an interesting alternative to clashing things together with higher and higher energy, at least a complimentary tool.

Another question has been experimentally adressed. What is the minimum number of atoms you need for a chiral system? It is one. Isolated Bi atoms have been shown to rotate the plane of polarised light.

I am sorry I am not quoting literature. The information comes from a lecture by Robert Berger, Stephen Hawking's "Brief History of Time", and Wikipedia.

[1] Two other fundamental symmetries are charge and time. Conversely to separate parts, it is generally assumed that CPT symmetry is conserved in the world, i.e. that things would be the same if you simultaneously changed charge, parity and the direction of time.

[2] Yes, I am still only treating physics at the popular science level. To be honest, it was even the Illustrated Brief History of Time ...

fold it

2008-09-07T19:45:00.006+01:00

Fold it is a pretty cool program I found at Lightnir's Blog is. If you are reading this blog it is probably because you are nerdy enough to look at 3D protein structures. If you are even nerdy enough that your idea of a computer game is playing with protein structures, you should check it out. If you do, look for FeLiXe tearing up the high score lists...

This what it looks like. Notice the amazing score of 9577 for this Calcium Ion binding protein. (I am leaving a little head start to lightnir ...)

The producers of the game mention that the computer game strategy may even be a fruitful way of tackling protein folding, which is basically the attempt to find the global free energy minimum structure of a protein. For every geometry you can compute the energy and forces. You can let the forces pull the geometry toward a local energy minimum. But the energy will usually be a very complex function of the geometry with many local minima. And this local minimum may be far away from the global minimum.

One strategy to go on is dynamics. Simulate protein motion at some finite temperature and hope that it will eventually overcome the barrier, leave the minimum, and go to a lower energy minimum.

The second strategy is Monte Carlo. Trying to systematically improve the structure by applying random changes.

The third strategy is making a game out of it and telling people it's fun.

Another interesting trivia about protein folding: Folding a random co-polymer is an NP-complete problem, meaning it has exponential scaling which makes it impossible to use for all but the smallest systems. The reason why in silico protein folding still works in some cases is first that we know the building blocks well and have a big knowledge base. The second reason is that proteins are folding in nature so we should be able to emulate this process and that there is a selection pressure for proteins to fold easily, so it works in nature.

The last thing to wonder about is the question between computer and human problem solving abilities. It will take something like 50 years until we can build something with the processing power of even an ant. Still we are using computers rather than ants for our everyday problems. Well it's because computers are more flexible. But it seems with anything that is remotely related to anything that evolution would select us to do we kick their asses.

Summer

2008-09-04T22:17:00.002+01:00

Alright let's go for a non science post. I just spent 8 weeks in the US. The bigger part at UIUC to do some research there. The idea of living on a campus is pretty cool. Just riding your bike to work. Stop by the gym during the day, tons of restaurants. We do not have that back at home. At home I always have to go through the big city with tons of traffic, tons of people,... But of course a big city is also cool.

On my last day I got to see all the ambitious incoming freshmen who still think that life makes sense and that school will take them somewhere. Maybe some even think that research is fun. Maybe because of the fact that it's fun to see how a little, easy sounding task takes up months until you finally have no idea about what you actually did. And when you finally kind of figure out what you are trying to say, putting that in writing takes even longer. Maybe I'll think about it differently when I finally have something published. Because then I can visualize those 10 or 20 people who will eventually superficially sift through the text ...
Well complaining about grad school is a nice tradition that should be upheld.

The Who's Who of Theoretical Chemistry ...

2008-08-14T03:39:00.004+01:00

... is basically this list. Well, of course a lot of great people have stayed back. Still it's amazing what names are on there. That feeds into my revolutionary spirit ...
I used to think science was about science, at least for basic research. Well even science is about squeezing every possible penny out of your competition.
Admittedly I am using Microsoft and admittedly I have used this program named after the 19th century German mathematician. I guess they are both ok, it's just sad how power makes you paranoid.

Energy

2008-08-09T23:01:00.011+01:00

It's been quite a while since I've written anything here. I don't know how exciting the topic is that I picked for coming back. But it's something to think about once. In physical chemistry (which is the common favorite of everyone doing chemistry) you often run across terms of the following form:

Now the question is how do we estimate that if we don't have a calculator at hand? And more important how do we do it with all the different units that are around?

We start out with a handy table that gives us the values of kT (at room temperature) in the different units that we are interested in. (Non-energy units are converted using Planck's constant and the speed of light in the usual way.)

298	K
2.478	kJ/mol
0.592	kcal/mol
0.026	eV
0.00094	Hartree
4.11E-21	J
207.1	cm^-1
6.21E12	Hz
48281	nm

Now we can estimate things because we know that:

and

For example we can say that a hydrogen stretch at 4000 cm^-1 will not be active at room temperature because it's energy is much higher than 207 cm^-1 which corresponds to kT. Or we can say that a reaction with a Gibbs free reaction energy below -25 kJ/mol will be almost quantitative.
But to be certain we need the following two formulae which will help to estimate the order of magnitude (ln is the natural logarithm, lg the decadic logarithm)[1]:

and

Now we can do pretty much everything we ever wanted to do.

In the first two examples we see that we are in the order of 10^-4 and the guess was correct.
If we stick to normal modes, we can for example be wondering to what extent a 600 cm^-1 normal mode is in its first excited state. Using the first formula we get:

That's below10% but well above 1% (the exact value is 5.5%). And we know that a mode at 600 cm^-1is still excited a little bit.

But of course you can also do it for more complex things. For example coalescence in NMR. For an isomerising molecule with activation barrier ΔG coalescence is reached if the following condition is fulfilled:

Δν is the frequency difference between the peaks, e.g. 0.5 ppm of 600MHz which is 300 Hz. kT/h is the frequency term in the Eyring equation, it actually corresponds to my value from above, 6.21E12 Hz, but I think for a different reason [2]. If you divide by it you get:

Apply the logarithm formula and put in the value of kT (or actually RT) and you get:

And that's about the actual value.

[1] The advantage is that you have to remember only to digits 2.3 and you almost get four digit accuracy compared to the exact value 2.302585...

[2] I would say that the value is chosen because it is a typical frequency of a skeletal mode but it is coincidence that this is close to kT.

Molecular dynamics movie

2008-06-23T19:29:00.006+01:00

A short instruction: how to make a movie out of a molecular dynamics trajectory with pymol. First you have to import the trajectory. The tricky thing is that pymol will not read multiple xyz's. But at Lightnir's I found that multiple sdf's work. So you start out converting your trajectory into sdf format (for example with openbabel). Then you open that in pymol, adjust the representation and save the movie frames. Either by something like:


set ray_trace_frames=1
set cache_frames=0
mclear
mpng mov/mov

or by doing that in the "Movie" menu.

Now you need a program that makes a movie out of the frames, for example videomach. And that's it.

Here's an example: the ground state zero point vibrations of ethene. Actually a classical simulation of this quantum phenomenon is neither physical nor exciting. But I don't have anything else to show right now.

Glucocorticoid receptor (3)

2008-06-17T22:12:00.003+01:00

Again the same molecule with a few more pictures. They were made with the "apple cutting" plugin from Lightnir. It cuts a plane through the protein. The orientation of this plane is not entered manually but chosen automatically by the viewing perspective which is a pretty cool idea in my opinion.

With pymol and this plugin you can quickly visualise a binding pocket and make pictures like this one ...

... or this one ...

... or this one.

Glucocorticoid receptor (2)

2008-06-08T13:26:00.005+01:00

I played some more with the receptor from last week. Here are some more images. The trick is setting multi-layer transparency in pymol (Settings - Transparency - Multi-Layer). Rendering takes quite long then, but it looks nice.

You can take a better lock at the ligand inside the binding pocket.

Or you can make this picture which looks kind of cool, I think.

This one shows the opening in the outside surface.

Two commands to mention are:

Changing the transparency of an object ("cha"):
set transparency, .1, cha
I can never quite remember the syntax for that.

Pymol has a nice way of parallelizing the rendering. The downside is that this will stop the whole system. If you want to keep working while rendering in a multicore machine, you can change the maximum threads:

set max_threads, 1

Glucocorticoid receptor

2008-06-01T11:06:00.005+01:00

Time for some more pictures. This (1P93 from the pdb) is the glucocorticoid receptor which would for example bind to cortisol. In this case it's a flourinated derivative of it. I did not read the story behind it but my point here is just to show a ligand in its binding pocket.

The whole protein looks like this, the ligand is shown in red.

If you draw the protein with its van der Waals surface, you notice that there is a small opening which is probably the place where the ligand entered the binding pocket.

Apparently the receptor closed again and the molecule is trapped in there.

If you eat away some of the outer amino acids, you can take a better look at the binding pocket.

Or at a semi transparent binding pocket:

Remove some more amino acids that are in the way and it'll look like this.

How do you create those picutres in pymol?

First you select one chain (the fact that there are four chains is just because that's how it crystallizes):
select cha, chain a
If you go for "preset pretty", you can see the ligand and select it. Call it lig. Then you can modify it, give it double bonds, show it as ball and stick model, ...

For the main protein you may want to select the different secondary structures:

select helix, ss h and cha2
select sheet, ss s and cha2
select loop, cha2 and not (ss h or ss s)

Then you can color those the way you want.
If you want a surface, it makes sense to duplicate that chain a and show the duplicate as a surface.
If you want to select the binding pocket, go for:
select pocket, (lig around 5)
Then you can erase some of the amino acids that are in the way.

Of course after doing all this, you have to play around with colors, perspectives and rendering methods.

NMR prediction

2008-05-23T16:17:00.004+01:00

Here is another pretty cool tool, prediction of ¹H NMR spectra.

This is the spectrum of the substance from last post as predicted there. It looks kind of like the real one and none of the shifts are off more than .5 ppm.

The original author's homepage has a better design but the second plugin is not working in my browser.

Generating structures

2008-05-22T17:09:00.007+01:00

I just heard of a pretty cool tool for generating structures.

Imagine the following problem. You have the structural formula of some molecule, let's say

And now you are wondering what the molecule actually looks like. If you have a "real" molecule model kit, you could use that. But you probably don't carry that around with you everywhere, it may not have enough atoms and it does not give you quantitative information and so on. The alternative is the Online SMILES Translator.

First you draw your molecule in a standard molecular editor (e.g. ChemSketch) and save it as a .mol file. Then you upload it to that tool and create a 3D mol file. Open the result in a 3D molecular editor: for example Pymol if you want nice pictures; Arguslab if you want some quantitative information and do computations.

The structure I just mentioned looks like this.

I did not really find any special information in this case. But for example it is interesting to notice that the part of the ring with the double bonds is planar.

The side remark is that there are of course more interesting things than just the global minimum (assuming that the tool even finds this minimum). For example other stable structures and the barriers to get a feeling of the dynamics. But I think it's still quite a bit of information for how fast it is and for how little input it needs.

Add on: At this site you can use the CORINA structure generating tool without the need of any other software. "Create molecule" opens a java applet editor, the result is shown in another applet, both in your webbrowser.

tRNA

2008-05-04T16:58:00.004+01:00

It's about time for some more pictures. I picked tRNA because of a course about modelling nucleic acids which I am doing these days.

The secondary structure looks like this, the clover leaf. The amino acid is bound at the top, the anticodon at the bottom. The lines represent the base pairing.

For the 3D structure you have to first consider that helices are formed with base pairing - twist around all those parts. In tRNA there is an important tertiary element, the connection between the left and right loop - fold the molecule. You will get something like this.

The open end is shown in read, the anticodon green, the right loop yellow, the left loop cyan.

Twisted around:

The open end:

The anticodon loop:

These three bases are probably the anticodon. They are GAA which correctly codes for Phe which this piece of tRNA is supposed to do.

If you want to move the molecule yourself you can find the link to a jmol applet from the pdb.

The place to be ...

2008-04-30T15:53:00.005+01:00

... this July is Vienna. Because then you can attend this pretty cool summer school:

If you are coming I invite you for a beer (unless you are coming either way ...)

By the way you can ask yourself who the two people on the picture are, and why.
And which one of them was on the Austrian 1000 Schilling note.

Some more exciting math

2008-04-28T17:41:00.010+01:00

I haven't written anything for a month and the thing that catches my mind is proving why the Hamiltonian operator is Hermitian, which is of course extremely exciting. What I liked about it, is the fact that the proof is not really more than the formula for partial integration.

A is Hermitian if

or equivalent

where <.|.> is a scalar product of a unitary vector space

As a warm up one can show that multiplying with any real constant λ is Hermitian. This is independent of the chosen scalar product and we can write (the last part is possible just because λ is real and therefore equal to its complex conjugate): [1]

Let's define the scalar product as:

Now comes the part that I mainly wanted to show, proving that 1/i d/dx is Hermitian with partial integration (please excuse that I am leaving out the lim's but I think it's easier to read this way)

partial integration leads to

The first expression vanishes because the wave functions have to be square integrable and therefore disappear at infinity. In the second expression we can draw the operator into the complex conjugation and the complex conjugate of i is of course -i.

And that's what we wanted to show. Actually a pretty amazing result if one weren't used to it.

For showing that the Hamiltonian operator is Hermitian you have to prove a few more things, e.g. that the square of a Hermitian operator is Hermitian.

You can do similar things with summation or real scalar multiplication.

Finally you can enjoy yourself with multiple integrals and find out that there is no problem with those either:

There is always only one variable involved in the operator and you can just add the other integrals on the outside. Of course this also holds for different x_i than x₁ because you can change the integration order.

[1] I am sorry I put a frownie face there.

Schur's Lemma (2)

2008-03-31T18:01:00.004+01:00

Here's the follow up post to Schur's Lemma, in case you were anxiously awaiting it.

The corollary says that if α=β, f has to be a multiple of the identity function (or unit matrix).

For the proof one looks at an eigenvalue λ_i. And subtracts the following on both sides:

This changes to:

Now we apply the Schur's Lemma from last time and find out that has to be either invertible or the 0 function. But it cannot be invertible because λ_i is an eigenvalue (actually the only eigenvalue). And we have what we wanted.

Quantum cryptography

2008-03-24T11:24:00.016+01:00

A pretty cool book I just read is Simon Singh's Code Book about cryptography. As my blog title suggests I'll mostly talk about the quantum part. Just a short summary.

Today's standard is the RSA cipher. Its speciality is that you need a different key for enciphering and deciphering. In other words: everyone can a encipher a message intended for you, but only you can read it. This is possible because exponentiating and subsequent modulo taking is in principle bijective (with an appropriate domain) but it is difficult to find the inverse function. At the heart of this method are two (large) prime numbers p and q. Their product N is made public and can be used for encrypting. For decrypting you need p and q. If p and q are large enough (and well chosen so that a few other mathematical tricks don't work), it is in principle impossible to deduce p and q out of N and only the person that initially multiplied p and q will be able to read the message.

Apparently quantum computers might be spefically suited for the task of factorising N. A regular computer would check numbers one at a time. The input for a quantum computer is a superposition of all the numbers to be checked. The computation can be carried out on this whole set. If done correctly the measurement at the end will lead the correct result with a high probability.

The way to escape from every decryption even quantum computers is quantum cryptography. It relies either on the uncertainty principle or on quantum entanglement. According to modern understanding of science it is in principle unbreakable. Unless of course some physicist eventually finds some of those hidden variables, how you can hop back and forth between parallel universes or some more realistic way of life beyond quantum theory.

P.S. There may be some quantum computers around already and people that know every dirty bit of secret about you.

Simplified molecular input line entry specification

2008-03-17T20:19:00.004+01:00

The whole name is almost as catchy as "SMILES". I used to think it was a strange way of representing molecules for computers. But it actually seems like a more straight forward way than IUPAC nomenclature. It's also shorter, and there's a possibility to have unique names. (It's more difficult to pronounce though.)

You can try out SMILES strings at this page it's kind of fun. How to do it is described on wikipedia for example.

Ethane is just CC.

Add double and triple bonds like this:
C#CC=C for butenyne.

Add a branch in parentheses:
CC(C)CCC for 2-Methyl-n-pentane

If you want a ring add a number after the two atoms to be joined together:
C1C(C)CCC1 for Methyl-cyclo-pentane

Add a pyridyl group to the C next to the methyl group (aromatic atoms are written in lower case, and you have to include a second ring closure)
C1(c2ncccc2)C(C)CCC1 for (2-Pyridyl-)-2-methyl-c-pentane

You can add an extra oxirane ring:
C1(c2ncccc2)C(C)CCC13OC3

You can mess with stereochemistry (using @ and @@)
C1(c2ncccc2)[C@@H](C)CC[C@]13OC3

If you still haven't had enough, you can add a double bond in E configuration to the pyridyl ring:
C1(c2nc(/C=C(Cl)\C)ccc2)[C@@H](C)CC[C@]13OC3

or Z configuration
C1(c2nc(/C=C(Cl)/C)ccc2)[C@@H](C)CC[C@]13OC3

Go SMILES!

Schur's Lemma

2008-03-09T09:15:00.013+01:00

Schur's Lemma is a nice piece of group theory apparently needed to prove the orghogonality relations which are in turn another nice piece of group theory. I did not find the whole proof anywhere, so I will show it here. The nice thing is: Once you put things into correct mathematical notation, you're almost done.

We have a group G whose elements we call R. And we may think of the elements as symmetry operations.

We have two representations α and β of this group. So what is a representation of a group? In most cases you may think of it as a vector space V or its basis, and let the group operations work on it. To be precise you have to define how the group operations work on this vector space. Then a representation becomes a mapping (a homomorphism) between the group G and the group of invertible linear transformations on V (instead of linear transformations you may also think of matrices).

Typically β=α but for this prove it can be different and may work on a different vector space W

For Schur's Lemma we also need that α is irreducible. That means every subspace U of V that is invariant with all the group operations contains either just the null vector or is the whole V.

Finally we have another linear transformation (or matrix) f

which has the following property

In other words the two mappings are the same for all the group operations.

Then Schur's Lemma says that f has to be either the 0 mapping or it has to be invertible. In the second case are the dimensions of the representations the same, dim(V)=dim(W):

If you have read this far, you may be wondering why that is.

Proof:

First we can assume that

without loss of generality because we could just relabel things if they aren't.

You have to consider the kernel of f. First we make the following observation

In other words:

But because α is irreducible, we know that ker(f), a subspace of V, has to look like:

In the second case obviously because everything is mapped onto the null vector. In the first case f is injective and it has to be bijective because W has at most the same dimension as V (it turns out that dim(V)=dim(W)).

That's what we wanted to show.

The applications of this are also kind of interesting. Maybe later.

Kilokayser

2008-02-29T16:58:00.001+01:00

I have a new favourite unit. It's the Kilokayser. For example a CH stretching vibration takes place at 3 kK.

Energy plot

2008-02-27T17:21:00.008+01:00

Let's get back to science. I have another little python script I would like to show. It's for drawing energy schemes with the help of Matplotlib. In my eyes it is faster this way than if you have to draw it by hand. The package can be downloaded from my homepage.

The output looks for example like this.
The input is a (rather simple) python script. If you want to use the package, the best idea is to modify the example scripts according to your needs. If you have a question you can ask me because I'd be happy if someone actually uses my package.

Here is another example that is a little bit more complex.

The God Delusion

2008-02-19T12:55:00.003+01:00

I finally got "The God Delusion" by Richard Dawkins. First I have to say it's a great read because of his writing style. If you don't think it'll disturb you, you should read it. If you think it will disturb you, you may want to rethink your beliefs unless you are happy as you are.

What really caught me though is Dawkins' fanaticism. It's not just about enlightening people and giving them new insight. Not about making religious groups more tolerant. No, Dawkins says that he wants an organised church of confessing atheists with a strong lobby (maybe with services of active non-worshipping). Something that almost made me collapse is in the foreword. Asked if he is a fundamentalist, Dawkins says no because the other side is simply wrong. Just think about that! I am sure he is a smart dude, but that seems way out of line.

The first thing we learned in philosophy at high school is the difference between reality and truth. Reality is what we perceive, truth is what is behind. To say that we perceive the true world is a religious statement. This is equally unlikely [1] as any other single relgious belief. If you don't have this belief, I don't think you can ever use the word "true" in scientific context. Science is just a tool to manipulate the environment but no more. To me most discussions come down to Socrates: you know nothing and if you admit it, you may be a small step ahead. Imagine you are plugged into a matrix (like in The Matrix) that allows for a proof of God. Then you would have this perfect proof of God but only in your perception [2].

In my eyes we don't need another intolerant religious organisation. Religion should be based on humility. Knowing that you know nothing [3]. You can't tell what is right. But why put materialism above everything else? As I know from Dawkins' "Selfish Gene", in the scientific view, it's all about genes trying to reproduce. Why does that have such a high value?

And since you don't know what is right, you may well choose a religion. Just don't hate and discriminate against everyone else. To me the problem with religion is not that people have wrong beliefs (we just don't know that). It's that groups are formed that hate and kill each other. But that's the same with patriotism [4], nationalism, tribalism, sports fan-ism and any other kind of group formation. I guess group formation comes from biology, just like cells join to form a body. It's a task for society to handle the interaction between groups just like it handles the interaction between individuals.

To get back to the topic, the short summary of what I read so far: I don't like how Dawkins promotes his own arbitrary school of thought but I like how he ridicules everything else.

And more important it seems that more focus should be used for applying evolutionary science on improving society rather than on trying to fight religious groups.

[1] It seems unscientific to me to talk about probabilities for things that are beyond our scientific reach like Dawkins does.
[2] You don't need evil robots for a scenario like this but I think it is a nice example.
[3] I may be a fundamentalist when it comes to "knowing that you know nothing" but the other side is simply wrong.
[4] I am glad I never had to pledge allegance to a flag (basically worshipping a coloured piece of cloth to help me kill people who are carrying a different coloured piece of cloth).

Rotation

2008-01-20T22:39:00.000+01:00

Angular momentum is something that took me a while to appreciate. But as one has to deal with it a lot (orbitals, rotating molecules, spectroscopy terms) it makes sense to take a closer look at it.

I think the problem is if one thinks too quantum mechanically. And then it's difficult to understand what's happening. Angular momentum can be represented by a vector that points in the axis of rotation, its length corresponds to the absolute value of the momentum. We can write:

Because of the uncertainty principle we cannot know all the components. In a spherical system, except of having all three bits of information, we have two: L² and L_z. This can be derived from the fact that the corresponding operators commute and must have a set of identical eigenfunctions. We cannot exactly determine L_x and L_y because their operators do not commute with L_z. In other words: We know how fast our particle rotates and we have partial information about the direction.

Actually quantum mechanics imposes a second restriction (but still: keep thinking about a rotating particle). L² and L_z are quantised, hence "quantum mechanics". We have

and

One thing you notice is that:

In other words there will always remain something for L_x and L_y and the total momentum cannot be clearly determined. But if you move to classical mechanics (large l), the difference will vanish.

In QM strong motion means a large derivative which in turn means nodal planes (if you consider real functions [1]). This works here, too: l is the number of nodal planes the wave function has on a given sphere. m is the number of nodal planes you pass when you move around the z-axis.

[1] Complex functions can oscillate without nodal planes. e^ix always has the absolute value of 1 but a non-zero derivative. It's a rotation in the complex plane.

Coulomb and exchange

2008-01-05T17:32:00.000+01:00

Let's keep going with what we did last time. It will get even more exciting than it already was.

Well I think it's nice to see what's behind some of those interactions and if I write it down here I'll find it again. But I admit that I usually read blogs with pictures rather than lots of weird signs.

The first expression was

Let us assume that 2 or more orbitals are different between the Slater determinants (pn-2). Then one of the terms in every product has to be 0 and the expression vanishes.

If p=n-1, then we have to claim the following:

This of course means that π=σ. There are now (n-1)! possible permutations that work with this. And the expression reduces to: (the n! cancels against the n! we started out with)

If p=n, we also need π=σ to have all the factors unequal to 0 (and equal to 1). The expression reduces to

The squared sign is of course 1. For every i there are (n-1)! peruations π with π(1)=i and u_i=v_i. Then we have

The second expression was

This time it turns out that one of the overlap integrals (and hence the whole product) equals zero if pn-3.

I'll only talk about the case p=n-2. If you want to have no overlap integral equal to 0, you first need

Aside from that you have to make sure that the mutually different orbitals are in the last integral. 1 and 2 have to be mapped onto n-1 and n.

As a next step you can acknowledge that π and σ have to either be the same or different by one transposition.

In the second case their signs are opposite. We get the following expression.

There are appropriately (n-2)!*2 permutations π that fulfill the requirements above and since the signs are opposite, this reduces to

The derivations seem to be similar for p=n-1 and p=n. As Levine tells me (who was in turn told by Parr) we get for p=n-1

and for p=n

(the n! is only there because I started with an extra n!)

The first integral is the Coulomb integral (with g₁₂=1/r₁₂). It corresponds to the Coulomb repulsion of two blurred electron clouds corresponding to the MOs. The second integral, called exchange integral, comes in because of the Slater determinant and is a consequence of the Pauli principle.

The next question is using the spin functions from two posts ago and see if some of the integrals vanish. And compare this to what you expected.

Slater determinants

2008-01-03T12:54:00.001+01:00

I am having some time these days to catch up on the basics [1]. This could therefore be another dry mathematical post. On the other hand: you may have been wondering all your life what exchange interaction was. Then this is your chance, in case you are open minded when it comes to looking at a few summation signs and angle brackets.

Wave functions are typically built on spin molecular orbitals, u_j(i), i=(x_i, y_i, z_i, m_si) (i.e. "i" is an abbreviation for the 4 coordinates from last post). MOs are functions that involve one electron. To get a multielectron function we could put each electron into an MO and multiply the functions:

In these functions electrons are statistically independent, therefore also uncorrelated (which is a weaker statement).

The product would violate the Pauli principle. In order to comply with it, we introduce a Slater determinant.

This is the explicit form of a determinant. You form the sum of all possible products with electrons in different orbitals with an appropriate sign (-1 or 1). (S_n is the symmetric group of all permutations; definition of sgn).

In a Slater determinant, electrons are not stastically independent anymore. But "correlation" usually means going beyond the formation of a Slater determinant.

Let's evaluate a matrix element between two Slater determinants.

All the orbitals are taken out of an orthonormal system. We will notice that the determinant vanishes unless most of the orbitals are equals. Hence we arrange the orbital in a way that equal orbitals are moved up front. Let p orbitals be equal.

The operator A shall be a sum of similar 1- and 2-electron operators:

Now we can write down the equation:

Because of linearity, we can take out the summations

Because of the antisymmetry, the integrals for all the electrons are the same

Now if you look at it as a multiple integrals, you will notice that you can take out all the functions that are not affected by the operator. (<...> means the integral over the respective configurational space)

You notice that most of the orbitals have to be the same, in order for those products not to vanish. How to go on and where the spin comes in, next time.

[1] It kind of helped that I broke my hand on my first snowboarding day this year ... There's still going to be february.

Spin

2008-01-01T16:11:00.000+01:00

It took me three books (Atkins, Kutzelnigg, Levine) to finally understand spin. Not its physics (I don't like physics) but what it does in quantum mechanics.

The problem is that you only have spin if you include relativity into QM (and no-one wants to do that [1]). Without relativity, spin has to be introduced as an axiom.

I think a simple way to understand it, is the following:

A one particle wave function θ(x,y,z) without spin depends on the three spatial coordinates. It's a function

The spin-assumption is that the function depends on a fourth coordinate, called m_s, that can take on the values 1/2 and -1/2 (or "up" and "down"). In other words Ψ(x,y,z,m_s) is a function

Integration over the whole configurational space is now understood as three integrations and one summation:

Of course you can go on with this and define an N particle wavefunction:

where

is the (differential) probability of meeting particle 1 at (x₁, y₁, z₁) with spin m_s1 while particle 2 is at (x₂, y₂, z₂) with spin m_s2, and so on.

What does a wave function with spin look like?

Since the typical approximation for the Hamiltonian does not include spin, we can assume the wave function to be separable:

where θ is the spatial function from above, and σ is a function whose domain consists of only two elements.

Apparently if we have two linear independent functions of that kind, we can make any other such function as a linear combination of those two. We choose two such functions called α and β with the additional requirement that they form an orthonormal basis (of the spin-function-space):

These conditions can be satisfied by the following (the 4 values have to form a unitary matrix):

Now θ(x,y,z)α(m_s) is called a spin orbital with α-spin. It vanishes unless the spin is 1/2 (if α(m_s) is chosen the way shown). With the α(m_s) and β(m_s) chosen orthonormally you can easily evaluate integrals involving spin orbitals (easy as far as the spin-part is concerned). But I don't feel like showing that today.

[1] Except of course for Dirac.

More mass = longer

2007-12-10T21:14:00.000+01:00

I hope you don't have ambiguous thoughts with this title. But that's the only title I could think about for this topic. And the topic is an excellent piece of math that will draw the crowds [1].

The Schrödinger equation, stripped of constants, looks something like this:

Now we switch to mass-weighted coordinates:

And we change the wave function to mass-weighted coordinates:

Now we don't need that tedious mass in front of the derivative because:

(Usually I would complain if someone shows a prove like this. But I think it's easier to see this way. And I also did it the "clean" chain rule way.)

Using this we get Schrödingers equation without the mass:

That means: the product of mass^1/2 and length is more natural than just length.
It also means: if an elephant walks for one meter it's the same thing as if a mouse walks 1000m (assuming that the squareroot of an elephant's mass is 1000 as much as the squareroot of a mouse's mass).

The idea is from this article. It turns out that the approach is handy if you want to propagate wave packets over a grid you obtained from linear interpolation. And who doesn't want to do that?

[1] Just like my diploma thesis never ceases to catch the interest of everyone I talk to.

Erythromycin

2007-11-02T23:00:00.000+01:00

I can see the attraction of creating weird looking molecules with lots of stereo centers. But, well, drawing them is enough.

This is Erythromycin's structure formula, as it is found on Wikipedia. It's a macrolide antibiotic used for people with an allergy to penicillin.

This is the 3D structure:

This is the 3D structure from a different perspective:

Nothing new, I know but I get most of my hits from people googling for such graphics. So here you go.

What is new though is this little brain teaser. You'll find a jmol applet next to the structural formula. The applet contains 8 models you can switch in between with your right mouse button. One of them is the correct structure the other ones are diastereomeres. Maybe I am the only one who thinks this is fun. But if you do try it out let me know if you found the correct structure.

And why is there always a need for new antibiotics? Because little germs are succesfully beating big pharma over and over again. My alternative tree hugging side might like that but of course my rational side would rather stay healthy.

Matrices

2007-10-07T12:19:00.000+01:00

Another little bit of math related stuff for all you people who like to think about it. It's from here and it took me forever to find out what they were talking about.

To start with (and to use this nice LATEX tool) what do we have:

A time dependent vector:

And we define the expected value operator <> which means averaging over time.

The question is what is represented by the following two expressions.

No, in this case the T's cannot just be ignored as I would usually do.

The first one is a matrix like

is a matrix.

The covariance matrix, in other words:

The second expression is a number, the sum of all the variances or the trace of the covariance matrix (which stays invariant with a similarity transformation).

As far as implementation goes, I can only stress how nice numpy is. All I had to do was the parsing. Numpy quietly converts my 66x19553 matrix into its covariance matrix and diagonalises it.

Multi reference analysis

2007-10-01T17:01:00.000+01:00

I have to summarize one more of those lectures from the conference. That in spite of the fact that readers aren't giving me very much incentive to do so ... Anyway, I'll talk about multireference analysis (MRA), again nothing very chemical but pretty nice and I want to review it.

The talk was from Janos Pipek. I don't find a reference but I'll write a short story about it.

Every starting chemistry student is shocked after hearing that in spite of having this nice variational principle everything is only based on AOs and those AOs are all based on Gaussians. The way out is wavelet MRA, for example with a Daubechie function:

It has many nice properties:

It becomes zero after a given distance.

Those functions next to each other are orthogonal.

It is easy to refine the basis if better results are needed for some areas. The large functions can be accessed as linear combinations of refined functions.

It seems that integrals can be easily evaluated.

That makes it a nice candidate for a basis function. They use it for jpeg compression and I am sure that they will soon use it in theoretical chemistry.

CESTC 2007

2007-09-28T19:46:00.000+01:00

CESTC 2007 is over. It was pretty nice, of course only for scientific reasons ... Anyway what I am going to talk about is only the scientific part. I want to cover a few lectures during the next days.

The first one has the advantage of being mostly mathematical and you don't need to know a lot about theoretical chemistry, but you need linear algebra. I always thought that in theoretical chemistry everything was hermitian or unitary and you didn't need to worry about singular value decomposition. But I guess I don't have an excuse anymore to skip that in my linear algebra script when I will do the exam.

It goes like this: When using configuration interaction singles we have to use a n_occ x n_virt (numbers of occupied and virtual orbitals) matrix C that represents the weights of all possible excitations.

where c_ij represents the weight of the component where excitation goes from the i^th occupied to the j^th virtual orbital.

This is not a very handy information especially when you try to visualise the orbitals.

But we know that for every real matrix orthogonal matrices U and V exist with the following property:

with

We only have to consider n=min(n_occ, n_virt) excitations any more. The orbitals they are taken from and given to, are the columns in U and V respectively. Typically only one singular value λ_i deviates much from zero. And the corresponding MOs are localised on the chromophores.

Give it up for linear algebra!

For more information read I. Mayer's article.

If you want some homework, you can think about how this works for a 4-dimensional CISD excitation tensor.

Still living

2007-09-21T16:47:00.000+01:00

What are readers trying to tell me that I get more hits than ever before now that I am not writing? Well, anyway I feel like supporting you with some more important bits of information. These days I am doing some summer research. I can't tell you what I am doing because you might steal the information, publish, and get famous. But I can tell you what it is like. I don't lead the easy life of synthetic chemistry. I don't get to set up a synthesis in the morning and wait around all day until it is done. Life in my field is much rougher [1].

The nice thing about computations is that no work that's older than 10 years is any competition. It's not the fault of the people from 10 years ago, they were alright. They just did not have the computers. I remember our first PC from that time was 20 times slower and 5 times as expensive as mine is now. For the same price computational power increased 100-fold. That's the difference between running five jobs over night or running a PC for a whole month. Thanks to anyone who is funding me through the computer game industry.

Next week I am going to the Central European Symposium for Theoretical Chemistry. That makes me feel pretty smart. I don't know how much I am going to understand. But I can spell the word symposium: symposium.

[1] A computational chemist sets up a computation in the morning and sits around all day until it is done.

Hückel package (2)

2007-09-05T18:26:00.002+01:00

I added a graphics function to the Hückel theory package that I have introduced in my last post. To some people visualising a graphical scheme is more meaningful than just numbers. If you are one of them: here you go.

This is the π energy scheme of benzene. Half of its orbitals are neatly filled with 6 electrons, 4*1 + 2 electrons or an uneven number of electron pairs, however you want to put it. On a second look you notice that the energy scheme is a hexagon, just like the molecule.

import hueckel
h_mol = hueckel.hueckel('***/benzene.mol',\
     calc_everything=True,   print_results=False)
h_mol.show_energy_scheme(can_width=300, can_height=300,\
                bg='white', arrow_length=.2)

If we take a carbon away (but no π-electrons) we have a cyclopentadienyl-anion. With 6 π-electrons, we have again the feeling of a filled subshell (or stability). A second look shows us that the energy scheme is a pentagon.

What does the energy scheme of cyclobutadiene look like? Yes - a square on its corner. But with 4 electrons we have the half-filled subshell which is only stable in transition metals, not here.

In fact it is true for every monocyclic conjugated hydrocarbon that the energy scheme looks like the corresponding regular polygon on its corner. Isn't that pretty cool?

But it kind of loses its magic once you put it in mathematical terms. For a monocyclic system of n C-atoms the energy of the k^th π-orbital is given as:

This corresponds to the polygon on its corner (in the scheme I set α=0 and β=1). Since the polygon is on its corner we need an odd number of electron pairs (4n+2) electrons to fill it properly. That's Hückels rule.

As far as I know there is nothing like Hückels rule for anything else than monocyclic systems. Both benzopyrene and biphenyl are totally happy with an even number of electron pairs.

By the way: all this was an ab initio calculation, I would say. No empirical parameters, just symmtery considerations and plausible simplifications. For more information you can read about the Hückel method at wikipedia or in a quantum chemistry book.

Finally you may ask yourself: Is a cyclopentadienium-cation stable? Let's try it out!

    h_mol = hueckel.hueckel('***/03-c-pentadienyl.mol',\
                  calc_everything=True, print_results=False)
  h_mol.set_total_charge(1)
  h_mol.show_energy_scheme(can_width=300, can_height=300,\
                               bg='white', arrow_length=.2)

Of course not.

The package can be found on my homepage under "Python scripts". When I was little I sometimes used to like writing German programs in Visual Basic. They can be found there, too. But there are two language barriers to overcome: German and Visual Basic.

Hückel package

2007-08-28T14:34:00.002+01:00

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

2007-08-22T16:14:00.000+01:00

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?

2007-08-07T20:45:00.000+01:00

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.

In-silico drug design

2007-07-20T15:21:00.001+01:00

This is the estrogen receptor. Women typically have a lot of the swimming around in their cytosol. Once estrogen binds, the complex moves to the DNA and transcription of the corresponding genes start.

In this picture there is only the ligand binding domain. They couldn't crystallise the whole protein because it was too flexible. Instead of showing different domains they chose to show the same domain three times. Alright: it's because of the crystallisation and XRD technique. The unit cell happens to have 3 unique molecules.

The small green things are the estrogens bound to their respective estrogen receptors. If you are trying to find a new drug that targets the receptor, these so called ligands and their binding sites are the most interesting part. The goal of in-silico drug design is screening a large library of molecules to find different molecules that also bind to a given receptor.

If you don't know what the structure of a protein looks like but you have the structure of similar proteins, you can go for homology modeling.

Once you have the structure of the binding site you try to dock ligands that will probably fit in your structure. First you have to find the right geometry of the ligand. After you have done that you have to score the docking, i.e. estimate the change in free energy of the docking process. This is rather difficult because there is no zero point energy to compare your docking pose with and because you have to include entropic effects.

Electrochemistry (2)

2007-07-18T19:43:00.000+01:00

Alright, electrochemistry (2) is coming up. Time to write about fuel cells. They are some pretty cool stuff in my opinion. And they are the most important research area in electrochemistry.

The net reaction taking place is pretty basic:
H₂ + 1/2 O₂ -> H₂O
You may probably have heard it before.

Other types oxidise hydrocarbons or alcohols with air oxygen. An important example is the direct methanol fuel cell that this lady at BBC likes to use to power her cell phone.

Using the reaction above for the production of electricity is not quite simple. The trick is that you need an electrolyte that conducts either H⁺ or O^2-.

If you want to conduct H⁺, the polymer electrolyte membrane fuel cell is your choice. The electrolyte is Nafion an ion exchange membrane based on teflon with sulfonic acid groups. Nafion only lets small cations pass (which is kind of opposite to typical ion exchange). With the help of platinum, hydrogen is stripped of its electrons at the anode. It travels through the Nafion membrane. It reacts with oxygen at the cathode and regains its electrons to form water. The PEM works at 80°-90°C.

You can also have the O^2- move around. One way to have that is ZrO₂ doped with Y₂O₃. Next to the yttriums there are oxygen vacancies. At 700°-1000°C its conduction is good enough for use in the solid oxide fuel cell. The solid oxide fuel cell is nothing for your cell phone, rather for a small power plant with good efficiency, especially combined with a Carnot process.

A third type is the molten carbonate fuel cell (300°C). Oxide ions combine with CO₂ to from carbonate. Carbonate travels through the fuel cell. CO₂ is released at the anode.

Alright that's it for now. I like my internship but coding in Python all day doesn't really make you want to sit in front of your computer after coming home.

Electrochemistry (1)

2007-07-13T19:33:00.000+01:00

I just had my electrochemistry exam. So I thought it's a good idea to give a quick summary over modern batteries and fuel cells. Both are systems that use spontaneous chemical reactions to generate electricity. Fuel cells work continuously, batteries discontinuously.

For every electrochemical process there are four components that have to be considered: anode, cathode, electrolyte, electron conductor. "Red cat" tells you that reduction occurs at the cathode.

Probably the most important secondary (i.e. rechargeable) cell ist the lithium ion battery. You find it in pretty much every cellphone or portable computer. If charged the anode consists of anionic carbon (LiC₆) and the cathode of lithium-cobalt(III/IV)-oxide (LiCo₂O₄). Both of those don't sound very stable. In other words we have a high voltage (about 3.5V). With high voltage and low density components the lithium ion battery has a very high energy density.

There is no way to use water in such a system (it decomposes at 1.2V). Instead ethylene-carbonate with LiPF₆ is taken (large PF₆^- means low lattice energy and therefore good solubilty even in an organic solvent).

The two half-cell reactions are:

anode: LiC₆ -> C₆ + Li⁺ + e^-
cathode: LiCo₂O₄ + Li⁺ + e^- -> Li₂Co₂O₄

Lithium goes through the electrolyte, the electron goes through the conductor (and your cell phone). They meet at the cathode and reduce Co(IV) to Co(III). The nice thing is that Li⁺ just intercalates back and forth without altering the electrodes. This "rocking chair" mechanism makes recharging possible. Compare this to the various types of lithium batteries containing lithium electrodes. It is not possible to make the electrode appear in its original shape through recharging. Dendrites will lead to short circuits.

Another way to get around dendrite growth is having the electrodes in their molten state. A typical example is the Zebra battery which works at 300°C and gives 2.6V. It has liquid sodium as an anode and NiCl₂ dispersed in liquid NaAlCl₄ as a cathode. The electrolyte is a solid sodium conducting ceramic, e.g. NaAl₁₁O₁₇. Zebra batteries are used in electric vehicles.

A more down to earth method is the nickel-methal-hydride battery (1.2V). We have an alloy hydride (MH) anode, a NiOOH cathode and a KOH electrolyte. Half cell reactions are (the proton is transported via the OH^-):

MH + OH^- -> M + H₂O + e^-
NiOOH + H₂O + e^- -> Ni(OH)₂ + OH^-

Fuel cells maybe tomorrow. Have a happy summer!

Retinal - Ultrafast is just about fast enough

2007-07-06T19:52:00.000+01:00

Ultrafast, barrierless processes have the remarkable feature that they occur in less than a picosecond and can therefore be modelled with molecular dynamics. The cis-trans isomerisation of photo-excited retinal is ultrafast. This is not primarily God's grace toward theoretical chemists (who like to model the process) but it ensures that the process is efficient and vision works. [1]

The trick with an ultrafast process is that after excitation, the reaction occurs (almost) without a barrier. The time is determined by the skeletal deformations, cf. a 333/cm vibration has a period of 100 fs.

To understand what is going on, it is sufficient to consider a smaller system than retinal, for example the pentadiene-iminium ion (the N instead of O comes from the protein). The system is strongly polarised toward the N. The first excited state (shwon in red) is polarized the opposite way.

After excitation (γ means "photon" [2]) the cis-doublebond becomes a single bond. Because of that it stretches out. Now there is an ecliptical single bond with weak π-conjugation. It spontaniously rotates.

The crucial step comes now: The excited state decreases in energy as the rotation progresses. The ground state's energy increases. At about 90° torsion angle their energies are the same. This is called a conical intersection. Now a radiationless transition occurs and we are in the ground state again. Ground state means double bond. And that means that a linear alignment is preferred. The molecule can either move back to the cis form or it keeps rotating the same direction and becomes trans. The ratio depends on the exact shape of the potential energy surface.

During that whole process we move down on potential energy surfaces. The geometry is relaxed to a minimum in the excited state. Then the molecule switches to the ground state where it has a local maximum. And the reaction keeps going until we reach one of the stable ground state structures. The energy needed to drive the process comes from the photon. But as soon as the photon is absorbed everything happens spontaneously.

Because of the short reaction time no competing reactions are to be expected. Fluorescence is with 10^-8 seconds 5 orders of magnitude slower. But also collisions that could take up energy aren't likely in that kind of a time frame. That means we can expect high quantum yield. It is 65% for retinal in rhodopsin and 25% for retinal in solution.

The information was mostly from this article by Garavelli (10.1007/s00214-005-0030-z).

[1] Sorry for making a religious reference in a scientific text. But where I come from the gap between religion and science isn't quite as big. Our religious people don't try to sabotage scientific theories for incomprehensible reasons and our scientists don't announce that they are atheists all the time. I don't know whose fault the situation in the US is. Maybe religious people should be nice and tolerant. Maybe scientists should be smart and on top of things. Maybe it's difficult for anyone to stop a fight once it has that scale and all you can do is defend your own interests.

[2] I don't like the hν for two reasons. First the energy is neither the only feature of a photon (it has momentum, spin, ...) nor is hν a unique identification for it. Second I don't like ν because it looks like a v and it does not abbreviate the word "frequency" like f does.

Rhodopsin

2007-07-04T15:05:00.000+01:00

The interesting thing about retinal is its photochemistry. I will talk about that soon. But this is another one of the protein-with-functional-group-picture posts. I still haven't quite realized how nice PyMOL is and I still enjoy making pictures with it. But of course also thanks to Röntgen and Bragg. And Edwards et al. who investigated the structure of this bovine rhodopsin.

This is the subunit with retinal (shown in red). Long stretched helices.

Retinal is bonded to a lysine residue (pink) to form an iminium ion. Or protonated Schiff base as they say if they are trying to confuse me. Retinal fits nicely into the helices. The π-system is almost planar only the double bond in the ring kind of sticks out. The 11-cis-bond (where I don't understand why it is called 11, probably a terpene nomenclature) is behind the green arc.

Carboxylic groups surrounding retinal are important for the charge distribution. You can select them in PyMOL with the following command (after making a selection "retinal" with retinal):

select SurrAc, (retinal around 10) and resn asp+glu

There is one at the iminium group ...

... and one at the β-ionone ring.

They give the molecule kind of a jungle feeling.

Quinine (for real)

2007-07-03T14:15:00.000+01:00

I wrote about quinine a while ago. I was told that I had the structure wrong. This is what it should look like.

The difference is the configuration at the stereogenic quinuclidine carbon.

An important difference to the "epi-quinidine" which I drew before is that a hydrogen bond is not possible. In the picture below it can be seen that the -OH and N are bent away from each other.

Rotation around the σ-bond that would lead to an H-bond is not possible for steric reasons.

This is apparent in the spacefill model.

I used the typical programs (I wrote about before) and the Open Babel package that I recently read about at Noel's blog. It's a convient way for converting chemical file types. It was very useful for passing structures between the different modelling and graphics programs.

Blog ethics

2007-06-22T20:35:00.000+01:00

My ankle is still not quite as good as it's supposed to be. Therefore I am still spending more time in front of my computer than I would like to. And with time distribution when studying, no not procrastinating, I got to visit more blogs than usually. Two important questions came to my mind.

The first one: Are blog writers with black blogs better people? I have to disappoint you: no not in general. Someone apparently suggested that at nOnoscience. The idea was that people should use Blackle instead of google. Showing a black screen should take up less energy than a white screen. But this is apparently only true for CRT screens, not for LCDs. So black bloggers: you are only better people if your readers have CRTs, not in general.

I guess a better approach would be just turning your screen brightness down.

The second question is even more important: Should I keep reading R-rated blogs? I found this important tool on a blog that is itself R-rated. It tells you what audiences your blog is appropriate for. The disturbing truth: most of the blogs I read are R-rated. It was no wonder with Ψ*Ψ/Excimer and Kyle. But even Mitch who I thought was responsible couldn't keep himself and his co-bloggers down. Even worse: I need parental guidance to read Albert's posts. And don't tell me words aren't as important as content. I am already disenchanted.

Hemoglobin

2007-06-21T21:15:00.000+01:00

I stumbled across hemoglobin while I was studying for toxicology. The body's oxygen carrier is a pretty cool molecule. Here are some pictures that I drew in PyMOL during little "studying breaks". They are made from the XRD structure from Protein Data Bank. Thinking about it, it seems pretty cool that we can actually visualise molecular structures and that we get the same results when we model them in a computer.

On this one you can see the whole tetramer. The four heme groups are shown in yellow.

Heme is the center-piece of the molecule (or as they say prosthetic group). The O₂ is right in the middle. Underneath is the Fe²⁺. It's bonded to the oxygen, the 4 nitrogens of the porphyrin ring, and the nitrogen of His⁸⁷ which is the covalent connection to the protein (hence prosthetic group).

A hydrogen of His⁵⁸ is very close to the O₂. Good for it that it is bent away. You could imagine that linear CO wouldn't be quite as comfortable in there. This is probably one of the reasons why the affinity of hemoglobin for CO is much lower than that of isolated heme.

This is heme with its whole subunit. You can see how the two histidines (the covalently bonded one and the steric hindrance) come from different helices.

Doing all this makes me wonder which I like better quantum chemistry or theoretical biochemistry. I'll get to take a glance into both of them this summer hopefully. Maybe I'll know more then.

jmol troubles

2007-06-15T13:03:00.000+01:00

Jmol is such a cool program that not using it in my blog is not an option. The problem is that it causes weird error messages. If you had problems accessing my blog the last couple of days, it's because of jmol. The weird thing is that it works alright if you click a link that leads you here. But if you type the address into the address line, it doesn't work. I know one thing: the problem wasn't with the jmolInitialize command. It doesn't seem to matter to call it several times on a page. The routine jmol.js was programmed in a way that it is just skipped if it's called a second time.

If you weren't able to access my blog, these are the posts that caused the trouble: Methane vibrations, p-Cresole vibrations. Take a look if you haven't been able to. And don't tell me they are not cool.

Cellobiose (2)

2007-06-14T21:15:00.001+01:00

I took a second quick look at cellobiose.

I was wondering if you could see the hydrogen bond (in the center) in the MO scheme. So I checked out the orbitals with a large coefficient at that hydrogen and found two with a strongly bonding interaction. MO 41:

MO 54:

The 53^rd MO doesn't fit in (it has even an antibonding interaction) but it looks kind of cool.

Ok, that's it I have to be studying these days. There are a few exams in the coming last two weeks of our semester that I want to take. A last outlook into different fields before it all comes down to theoretical chemistry next year. I finally have to say that thermal engineering is kind of cool with all its balances. As long as you don't actually have to build the machines.

Cellobiose

2007-06-12T19:30:00.000+01:00

Cellobiose is interesting because it is the building block of cellulose. As I understand it, the important structural feature is the hydrogen bond between the ring O and the 4' OH. It gives rigidity to cellulose strands.

In my GAMESS 3-21G calculation the hydrogen bond has a bonding order of .08 which is not really a lot. But this may just be a problem of the calculation. I don't know how to add extra functions just to the atoms I am interested in. And even the way I did it, it took 3 hours of geometry optimisation and it wasn't even quite stationary then.

This post is mostly to show a few nice pictures, though. I like making them and it seems that most of the google hits are going for them. Here's another one.

To give credit to all the nice free (or formerly free) programs I used: The structure was drawn, MM optimised and semi-empirically optimised in ArgusLab. Then I rearranged the structure with PyMOL's sculpting function to get the O and OH together. Structure optimisation in GAMESS, data extraction in ChemCraft, pictures again in PyMOL.