Molecular Modeling Basics

Negative activation energies in molecular modeling: diagnosis and cures

2012-01-22T11:52:00.001+01:00

This post is inspired by a recent discussion with Paolo in the comments section of this post.

So, you've found that the energy of your transition state (TS) is lower than your reactants (i.e. you have a negative activation energy). Don't panic, instead pretend you're Dr House and go through it methodically:

Did you actually find the right TS?
A TS is a completely optimized geometry (no constraints, "zero" gradient) with one has one and only one imaginary frequency. Yes? OK, how big is the imaginary frequency and does the normal mode look like what you would expect? Sometimes minima can have small (< ca 100 cm-1) imaginary frequencies due to numerical "noise", or the TS search algorithm will find a TS for, say, methyl rotation. No, looks OK? (Consider computing an IRC do be really sure, just in case.) Let's move on.

How much lower?
Some processes have very low (< 1 kcal/mol) barriers. Numerical "noise" can cause these barriers to come out slightly negative (-0.1 to -1.0 kcal/mol) instead. If so, consider your reaction barrier-less for all intends and purposes. No, much more negative? Let's move on.

What energy are we talking about here?
If the activation energy is calculated using single point energies (for example B3LYP/6-31G(d)//RHF/3-21G) then the B3LYP/6-31G(d) PES may not have a barrier or the B3LYP/6-31G(d) TS geometry looks very different from the RHF/3-21G TS geometry.

The easiest way to check for this particular case is to geometry optimize your reactant at the B3LYP/6-31G(d) level of theory. If the optimization results in the product geometry, then there is (very likely) no barrier to the reaction. If the optimization gives you a reactant geometry, then one explanation is that 3-21G is not a reliable method for finding the TS and the solution is to use B3LYP/6-31G(d) in your TS search. There are other explanations, so read on.

If the activation energy is calculated using zero point energy or free energy corrections and is negative (this is not uncommon for reactions whose PES barriers are less than 2-3 kcal/mol) then you can consider the reaction barrier-less. If you want to be really sure, you can add these corrections along the entire IRC. However, there are also other explanations so read on.

Is your reaction unimolecular?

R -> TS -> P

OK, so you have found a reasonable-looking TS using some method and the electronic energy barrier computed using that method is quite negative. If your reaction is unimolecular then the most likely explanation is that you have not found the lowest energy conformation of your reactant. If can't find a lower energy structure for your reaction, compute an IRC. An IRC follows the minimum energy path downhill to your reactant and product, so you will find structures with a lower energy than your TS.

Is your reaction bimolecular?

R1 + R2 -> R1/R2 -> TS -> P

OK, so you have found a reasonable-looking TS using some method and the electronic energy barrier computed using that method is quite negative. The first question is how did you compute the barrier?

If you computed the barrier as the energy difference between the TS and the bound complex of your two reactants (R1/R2) then the most likely explanation is that you have not found the lowest energy conformation of this complex (see the section on unimolecular reactions).

If you computed the barrier as the energy difference between the TS and the sum of the energies and your two isolated reactants [E(R1)+E(R2)] then there is not necessarily anything wrong with your calculations. Read on.

How should I compute the activation energy of a bimolecular reaction?
If you are trying to model a reaction in solution your activation energy is

E_a = E(TS) - E(R1/R2)

This is because the energy liberated by binding is lost to the environment before it can be used to react. This activation energy should not be negative (see Is your reaction unimolecular?).

If you are trying to model a reaction in the gas phase your activation energy is

E_a = E(TS) - [E(R1) + E(R2)]

If this energy is negative then you would say that this is a barrier-less reaction and the reaction rate is proportional to the collision frequency. This is because the energy cannot escape the complex, and will eventually be used to overcome the barrier. However, this will be pressure dependent (higher pressure, worse assumption).

Animations of Unseeable Biology

2012-01-14T09:58:00.002+01:00

More animations by Drew Berry

Related blog post
The computational Microscope

The ideal app for an ideal gas

2012-01-08T11:33:00.000+01:00

I recently came across a very cool app called Atoms in Motion ($2.99). With this app you can perform molecular dynamics simulations of He, Ne, Ar, and Kr (and mixtures thereof) at different temperatures, volumes, and number of atoms.

It is effortless to set up a simulation and change the conditions: swipe to change T (or shake the iPad!) and pinch to change V. You can also flick an atom and, frankly, the most fun is to build a few clusters at 1 K and then move an atom carefully in position and flick it at a cluster! Something I have dubbed "atomic billiards".

Teaching
Since not every student has an iPad (yet), the main teaching-use will probably be in lecture, projection the simulation on using a VGA adapter (hint: simulations make for great peer instruction questions). Because it is based on real simulations, and not pre-recorded movies, the app is very versatile and can be used to demonstrate very simple things like solids, liquids, and gases or an atomic view of temperature or relatively complex things like defects in solids and deviations from ideal gas behavior (using the plotting tool).

However, the ideal use of this app would really be self-guided discovery. Most features are very intuitive (the plotting tool takes some practice) and I'd love to see what a class of 10-year olds iPad users would get out of this. Side note: it would really have been great if the iPad could vibrate and the pressure could be linked to the vibration.

Documentation
Another strength of this app is the quality of documentation and extra information available in the app and on the top bar of the website. I'm particularly impressed by the animated figure on "how big is an atom", which reminded me of the orders-of-ten approach. I would be lovely to see a similar explanation of the nano-second.

Disclaimer: I was made aware of this app by the developer, but I bought my own copy, and was not asked to write a review on the blog.

Bond lengths, Avogadro, Google Docs, and gamification

2011-12-10T11:00:00.000+01:00

Recently a colleague of mine, Henrik, lamented that at a recent oral exam many first year students had no clue what a typical bond length in a molecule is. I am not a big fan of memorizing something that can be easily Googled, but we're talking order of magnitude here, not the last decimal place.

This is a place where molecular modeling tools such as Avogadro can help, but how to do this exactly? Knowing an approximate bond length by heart comes from looking at many, many molecular structures, and imagining a worksheet with dozens of fill-in-the-blank bond lengths already puts me to sleep.

But how about something like: use Avogadro to build a molecule with the longest CC single bond possible. This sounds a bit more fun and I bet a few strategies are already popping into your mind.

Based on my recent excellent experience with Google Docs, I would then have student collect answers on a common Google spreadsheet. Students could either collect their work on a single sheet or one sheet per student/team. This allows students to learn different design principles from each other and adds an element of competition (gamification) since you can now compare results.

Would this actually work? I challenge you to find out!
By work I mean: engage students enough to build a large number of molecules and think about bond lengths? I don't know, but how about testing it out here? I have started a spreadsheet with a few entries. If you can do better, feel free to add your entries.

Some ground rules:
1. The bond length must come from the lowest energy conformation
2. To beat a current entry the CC bond length must be longer by 0.005 Å
3. You must use the MMFFs force field (which will limit the number of elements at your disposal).
4. For now I'll leave the maximum number of atoms open, but the fewer atoms you use, the more impressed I will be.

You might have the students figure these issues out by themselves. Explaining why they are important could be very teachable moments.

Practical tips
1. As the molecules get more complicated, naming them becomes difficult. In that case I suggest using SMILES, which you can generate with this GUI. It would be nice if Avogadro has a "Copy SMILES" option.
2. You can build your molecule using SMILES.

Using Google Docs spreadsheets when teaching molecular modeling

2011-12-06T16:04:00.003+01:00

Every year I teach DFT in the Computational Chemistry course here at the University of Copenhagen, where I give a 2 hour lecture one day, and give them an assignment to complete another day in the computer lab across the hall from my office.

I typically start them off around 10 am and they typically finish around 2 or 3 pm, when I then lead them through a discussion of the data. In past years this has been pretty chaotic: computed energies are missing or wrong, and relative energies had to be converted to kcal/mol on the fly.

This year, inspired by Luca De Vico, I set up a Google spreadsheet where they could collect their data (they work in pairs) on separate sheets. This worked great:

The students where able to see each others data, and they quickly discovered when they made mistakes (which I then helped them fix).

A "best practice" with respect to organization and presentation of the data was quickly adopted.

Discussion of the data was made much easier since I could just show the spreadsheet on the projector.

Finally, I could follow their progress from my computer and figure out when most people had finished, so we could start discussion.

Highly recommended!

Molecular Modeling Basics Reviewed in ChemPhysChem

2011-12-04T09:15:00.000+01:00

A very nice review of the book (and mention of the blog) by Ralf Tonner in ChemPhysChem. Much appreciated.

Reinventing Discovery: Practical steps toward open science

2011-12-01T08:16:00.001+01:00

What can you do if you're a scientist?

I recently finished reading Michael Nielsen's book Reinventing Discovery. This is not another review of the book (there are many: Google it). I'll just say it is a well written book on a very important topic, so go buy it now.

Towards the end of the book there is a section entitled "Practical steps toward open science" and one subsection is about what you can do if you're a scientist. There are some good suggestions, some of which I'd like to expand on here, and I'd like to add some new suggestions as well. The suggestions are roughly ordered in increasing effort (though not necessarily impact!), together with my own examples (where applicable).

Read the rest of the post over at Proteins and Wave Functions

"Missing" MOPAC parameters

2011-12-01T08:14:00.001+01:00

A long, long time ago in a land far, far away I implemented MNDO, AM1, and PM3 in GAMESS. This was done by taking chunks of code from MOPAC that contained contained the parameters, integral code, and Fock matrix builder. In doing so, I never noticed that the parameter file contained more parameters than where published in the papers describing the method. This conundrum came back to haunt us recently as we're trying to implement PM6.

Read the rest of the post over at Proteins and Wave Functions

Meanwhile, in another part of the blogosphere

2011-10-30T10:38:00.000+01:00

From Henry Rzepa's blog
Computers 1967 - 1985
Computers 1985 - 1989
Computers 1990 - 1994
Moore's Law and molecules
Related: Steve Jobs and chemistry: a personal recollection

Video lectures by Luca De Vico on Proteins and Wave Functions
Basis sets
Optimization techniques - part 1
Optimization techniques - part 2

Peer instruction questions: internal energy

2011-10-09T12:54:00.000+02:00

Peer instruction questions: internal energy

View more presentations from molmodbasics

These slides show questions I used when teaching internal energy using peer instruction. The slides are in Danish, but I hope you get the idea and there is always Google translate. Any questions, just leave a comment.

The slide use a Molecular Workbench simulation you can access here (or download here once you have MW installed), and is further described in this blog post.

I start and carefully explain the simulation on the left, then ask what we will see when we start the simulation on the right. Then vote, and explain answer. Then finally run the simulation on the right.

Related blog posts
See all posts related to peer instruction here.
Internal energy and molecular motion
Simulations in teaching physical chemistry: thermodynamics and statistical mechanics

Chemistry of the future: 3D and augmented reality

2011-09-26T21:13:00.001+02:00

Way cool! 'Nough said.

I have a cold: a molecular perspective

2011-09-15T18:51:00.000+02:00

Correction to Eq 2.45

2011-09-04T12:16:00.002+02:00

Astute reader James McNeely found an error in Eq 2.45. Here is the correct version. Thanks James!

I have started a separate page listing known errors here. If you find mistakes in the book please leave a comment on that page (or send me an email).

Rotating molecules in Powerpoint - part 2

2011-08-14T12:27:00.001+02:00

Animation in Powerpoint
One of the most accessed posts on this blog is Rotating molecules in Powerpoint. I wrote the post more than a year ago, and had I written the post today it would be different thanks to a great free program called Screencast-o-matic.

The screencast above shows how to use Screencast-o-matic to make short movies of rotating molecules and other molecular animation and include them in Powerpoint presentations.

Jmol in Powerpoint
A lot of people end up on this blog by searching for "Jmol in in Powerpoint" or some similar term. As far as I know it is not possible to embed Jmol in Powerpoint slides. It is possible in Beamer, as discussed in this post by Janus Eriksen.

Using Screencast-o-matic you of course record any sequence of Jmol animations, but it will not be interactive. Personally, I make a web page with my Jmol model, with buttons to control what I want to do and simply switch between Powerpoint and a browser. I give a short example using this Jmol model in the screencast below.

Blurring the boundary between linear scaling QM, QM/MM and polarizable force fields Part 2

2011-07-24T13:39:00.005+02:00

My talk at WATOC 2011. Slides can be found here

Here's how I made the video: I recorded my talk using the Voice Memo app on my iPhone. Then I replayed the talk on my Mac using iTunes as I went through my slides on Powerpoint. I recorded the screencast + audio using Screenflow.

Thanks to Anders Christensen for recording the video.

Summarizing a paper using Prezi and Screencast-o-matic

2011-07-21T15:06:00.001+02:00

In a previous post I summarized a paper using two programs: Prezi and Screencast-o-matic. Both are free and easy to use. The screencast above shows how I did it. Anyone can do this.

I did this on a mac, so I used the earphones/microphone that came with my iPhone. The sound is not the greatest, but good enough I think.

The screencast is 9 minutes long, which might be too long. It's hard to strike a balance between detail and the big picture. Hopefully, I will improve with time. Feedback is very welcome.

Update: be sure to check out the comments below

New paper: Ring current effects in proteins

2011-07-21T08:03:00.000+02:00

Definitive Benchmark Study of Ring Current Effects on Amide Proton Chemical Shifts

Anders S. Christensen*, Stephan P. A. Sauer, and Jan H. Jensen*

Department of Chemistry, University of Copenhagen, Universitetsparken 5, DK-2100 Copenhagen Ø, Denmark

J. Chem. Theory Comput., 2011, 7 (7), pp 2078–2084

DOI: 10.1021/ct2002607

Abstract (the paper is summarized in the video at the end of the post)
The ring current effect on chemical shifts of amide protons (Δδ_RC) is computed at the B3LYP/6-311++G(d,p)//B3LYP/aug-cc-pVTZ level of theory for 932 geometries of dimers of N-methylacetamide and aromatic amino acid side chains extracted from 21 different proteins. These Δδ_RC values are scaled by 1.074, based on MP2/cc-pVQZ//B3LYP/aug-cc-pVTZ chemical shift calculations on four representative formamide/benzene dimers, and are judged to be accurate to within 0.1 ppm based on CCSD(T)/CBS//B3LYP/aug-cc-pVTZ calculations on formamide. The 932 scaled Δδ_RC values are used to benchmark three empirical ring current models, including the Haigh–Mallion model used in the SPARTA, SHIFTX, and SHIFTS chemical shift prediction codes. Though the RMSDs for these three models are below 0.1 ppm, deviations up to 0.29 ppm are found, but these can be decreased to below 0.1 ppm by changing a single parameter. The simple point-dipole model is found to perform just as well as the more complicated Haigh–Mallion and Johnson–Bovey models.

The paper Instructions on how I made the video will appear in a future post. In the mean time, enjoy!

Blurring the boundary between linear scaling QM, QM/MM and polarizable force fields

2011-07-15T12:30:00.001+02:00

Jan Jensen's presentation at WATOC 2011

View more presentations from molmodbasics

Thursday morning, room B

Update: a recording of the talk can be found here

Peer instruction: mixing

2011-06-19T11:46:00.001+02:00

Peer instruction questions: mixing

View more presentations from molmodbasics

These slides show questions I used when teaching mixing functions using peer instruction. The slides are in Danish, but I hope you get the idea and there is always Google translate. Any questions, just leave a comment.

Some comments about specific slides:
Slides 1-3 show results from a Molecular Workbench simulation, which you access here. It is a variation on the simulation I used to illustrate entropy. If you understand why the gas expands in that example, you also understand why the gasses mix.

Slides 4-10 show results from a Molecular Workbench simulation, which you access here.

Slides 4-5: after two votes there was no clear consensus, but most students likes A "no interactions between molecules". The second vote came at the end of the first of three back-to-back lectures, so we had a third vote after the break. There really was intense discussion of this during the break, when most finally settled on B "equal interactions between molecules". I think A was popular because "ideal solution" conjures up an analogy to "ideal gas", but how can you have a liquid if there are no intermolecular attractions?

Slide 7-8: here is alternated between the simulation and the question a few times. A better approach would have been to include the question on the web site with the simulation.

Related blog posts
See all posts related to peer instruction here.
Illustrating mixing
Simulations in teaching physical chemistry: thermodynamics and statistical mechanics

Peer instruction: radial distribution functions

2011-05-19T13:26:00.000+02:00

Peer instruction questions for radial distribution functions

View more presentations from molmodbasics

These slides show questions I used when teaching radial distribution functions using peer instruction. The slides are in Danish, but I hope you get the idea and there is always Google translate. Any questions, just leave a comment.

Some comments about specific slides:
Slide 8: The hint is given after the first vote.

Slides 11-18 show results from two Molecular Workbench exercises, which you can download here and here, once you have installed Molecular Workbench on your computer.

Slide 11: First I run the solid simulation, then pose the question and have a couple of votes, then click on the "show pair correlation function" in the MW simulation. Note that you have to run for at least 100,000 fs to get good statistics (i.e. relatively smooth curves). Then I explain the answer (slides 12 and 13)

Slide 14: Same procedure as slide 11, but for the liquid.

Slide 16: No vote, since the answer is pretty obvious, but I do ask where the peak at r = ~1.5 Å comes from. Of course it comes from the attractive part of the Lennard-Jones potential, and you can clearly see some particles sticking together in the gas simulation. To check, I change the depth of the well from 0.1 eV to 0.001 eV (simply double-click on any of the particles, and you will see what to do), re-run the simulation, and show the radial distribution function (summarized in slide 17).

See all posts related to peer instruction here.

Peer instruction: entropy

2011-05-18T14:43:00.000+02:00

Peer instructions questions on entropy

View more presentations from molmodbasics

These slides show questions I used when teaching entropy using peer instruction. The slides are in Danish, but I hope you get the idea and there is always Google translate. Any questions, just leave a comment.

The slides refer to two Molecular Workbench simulations, which you can access here and here, and read more about here and here.

Specific comments to the slides
Slides 1-5: Here I start the simulation without removing the separator. Then pose the question, two rounds of votes (usually not needed for the first simulation with 100 particles, as more than 80 % submit the correct answer), then remove the separator and explain.

Slides 10-15: Here I run the first simulation (small volume, low temperature) and carefully explain what the recorded times mean and how they correlate with probability. Then I pose the first question (slide 10), two rounds of votes, then run the double volume simulation, then summarize the right answer (slide 11), and explain it (slide 12). Same process for doubling the temperature.

See all posts related to peer instruction here.

Related blog posts
Illustrating entropy
Entropy, volume, and temperature
Simulations in teaching physical chemistry: thermodynamics and statistical mechanics

Peer instruction: force fields and energy minimization

2011-05-16T09:46:00.001+02:00

Peer instruction questions on force fields and energy minimization

View more presentations from molmodbasics.

These slides show questions I used when teaching force fields (slides 1-4) and energy minimization (slides 5-8) using peer instruction. The slides are in Danish, but I hope you get the idea and there is always Google translate. Many of the equations are discussed in Sections 1.5 and 2.1 of Molecular Modeling Basics.

I also use the Avogadro program when teaching force fields. See for example this post.

I have made a Jmol model to illustrate energy minimization, as discussed in this post.

See all posts related to peer instruction here.

Molecular modeling and peer instruction

2011-05-14T12:29:00.001+02:00

Last Tuesday I had my first go at using peer instruction. Peer instruction is a fancy term for a relatively simple process:

1. Put up a multiple choice question
2. Ask the students to come up with an answer without talking to each other. Give them 2-3 minutes.
3. Have a vote (more on how below)
4. If 80 % or more get the right answer, show the right answer, give a short explanation, and move on.
5. If less than 80 % get the right answer, ask the students to discuss it with their neighbor. Give them 3-5 minutes.
6. Have another vote.
7. Give them the right answer and explain it.

Voting using Polleverywhere.com
The screencast (best viewed in full-screen mode) shows how I use the website Polleverywhere.com for the voting. It requires the students to have access to a browser, so I ask them to bring a smartphone or a laptop to class. It is also possible to submit votes by SMS, but Polleverywhere only provides a British phone number for Europe, so this gets pricey in Denmark.

The free version of Polleverywhere.com can register up to 30 votes. For larger classes your department needs to buy a six or 12 month subscription, the latter is about $900. You just create an account on the site, create your poll, and you are ready to go.

All my questions have four options labelled A, B, C, and D, so I just create one poll per course. The example in the screencast is used in a course called KemiF1, so the codes for the options are kemif1A - kemif1D. You do need to create completely unique codes. The site assigns you numerical codes by default, so you can also just use those. The website with the question is here. Option D is "don't know". It is important to have that option so that everyone casts a vote every time.

I thought it would be cool for the students to see the votes coming in live, but they felt it biased their vote too much, so it is better to leave the question up as they vote, and show the results when the vote is done. You can monitor the progress on the vote on another device such as a smartphone or iPad.

Why use peer instruction?
Much has been written about peer instruction (just google it), but I find the two videos posted here especially informative. And think about this: if you lecture on organic chemistry and less than 80 % of your students understand what you mean when you draw a hexagon, is there really any point in going on to more complicated material?

Peer instruction and molecular modeling
The real challenge when using peer instruction is to come up with good questions and I think molecular modeling and visualization can contribute a lot here. For example, the spinning models in the example reinforces the fact that cyclohexane is not planar. Note that the students can go to the page and interact with each model as they mull their options.

I will show more examples of using molecular modeling and peer instruction in future posts.

Interactive chemistry ebooks: interactive figures

2011-05-09T08:11:00.000+02:00

Here is version 2 of my ebook mockup page (see here for a description of version 1). Actually there are two versions: one using Jmol and one for the iPad.

The Jmol version is how I would like the page to look. The static figures are replaced with interactive Jmol versions embedded right in the page, and the figure captions contain links (denoted by ">>") to larger version of the figures on a separate page, that also contain extra features such as animations (Figure 3.1) or overlays (Figure 3.2). An alternative design would have been to keep the static images, but insert links to the larger figures. I think I like the current look better. But for pages with more figures loading-time could become an issue with Jmol.

The iPad version is what is currently possible for the iPad with ChemDoodle Web Components. We are waiting for two things to happen: 1) WebGL in Mobile Safari, which will allow real 3D representations of the molecules, and 2) Implementation of contours and surface maps in ChemDoodle Web Components.

For the iPad version I generated the necessary ChemDoodle Web Component HTML code using the ChemDoodle software package, as I show in this screencast. The mol2 file I start with is generated with Avogadro.

ChemDoodle costs $60, and normally I restrict this blog to free software, but I am making an exception here because I think iChemlabs deserves to be supported since they made ChemDoodle Web Components open source. Furthermore, since I have access to the source code I could probably figure out the format for the HTML code, if I wasn't so lazy.

Related blog-posts:
Interactive chemistry ebooks: highlight and annotate
Interactive chemistry ebooks: let us start now!
ChemDoodling: now in 3D, but not (yet) on the iPad
ChemDoodling on the iPad and the future of interactive chemistry text books
iPad: even 3D molecules that can be viewed from any angle

Interactive chemistry ebooks: highlight and annotate

2011-05-08T11:07:00.001+02:00

In a previous post I argued that the best way to get started on making interactive chemistry ebooks is to do it on the web with straight HTML, rather than wait for the epub technology to catch up. But how would one highlight text and scribble notes in such a book?

One solution is to use Diigo, which is a free add-on to your browser as well as a free iPad app. The picture above show a screenshot of annotation I added to the web page mock-up I talked about earlier. Notice the "Web Highlighter" in the bookmarks bar, which, when clicked, gives you the blue tool bar you see. The highlights and annotations are stored in the cloud, so they can be seen and modified from any browser where you have installed Diigo. It also appears possible, though I haven't tried it myself, to share your notes with a group of people you define, which sounds like a very interesting study tool.

The computer version of the Diigo browser plug-in also has a "read it later" option, where the web page is saved and can be read in the Diigo app off-line. For some reason this option is not there when using the iPad browser. It also not possible to annotate in the Diigo app. But I wouldn't be surprised if both options would appear in future releases.

Related blog-posts:
Interactive chemistry ebooks: interactive figures
Interactive chemistry ebooks: let us start now!

Interactive chemistry ebooks: let us start now!

2011-05-05T08:30:00.004+02:00

The problem
Henry Rzepa writes thought-provoking blogposts and his latest, entitled What is the future of books?, is no exception. He describes the current state-of-the-art in science e-publishing and the very limited options and steep learning curves facing aspiring authors in this area.

Currently, there is only one publisher of interactive science e-texbooks: Inkling. As a result, Inkling has their hands full and is unlikely to be interested in anything but the biggest best sellers in this area.

This means we chemists have to go at it alone, but there a is a big practical hurtle: there a no tools available to us. The ePub 3 format that will allow interactivity will be released this month. But when will tools to create and read ePub 3 files become available?

My first, and so far only, attempt at creating an ePub document taught me that the e-editors and e-readers are completely geared towards fiction. Looking at Chapter 3 of Inklings Biology textbook on the iPad, which is available for free, is very educational. Their design looks nothing like the traditional book that current e-readers such as Stanza and iBooks try to emulate with page swipes and similar layout considerations. So we are not really waiting for these apps to become ePub 3 compatible; we are waiting for a science e-book reading app. Who will write this and when?

One solution
But do we need all this stuff? Broadly speaking, ePub is HTML in a zip file and ePub readers are hacked web browsers, and both are aimed at letting you read HTML off-line. If we assume internet access, the need for ePub is mostly gone and this allows us to get started now by doing what we already know how to do: make interactive web pages. Even if you think off-line access is essential, I hope you'll agree that "web-books" allow to test some design ideas for science e-books in ePub 3 format or book apps.

To kick things off, I have made a web page mock-up of how the section on electron density from my book Molecular Modeling Basics might look in the hands of the Inkling people, using an iPad-friendly template kindly made available by Matthew James Taylor.

This is version 1, quite bare bones, and with a few blemishes. For example, the top bar moves on the iPad (and on Chrome) because Mobile Safari does not respect position:fixed, and the Table of Content menu at the top doesn't work. I would very much appreciate tips and suggestions from HTML experts on how to improve such things. You are very welcome to copy the HTML and the associated css, and experiment with this yourself. (I am quite confident that a skilled HTML programmer could make a web page that is essentially indistinguishable from a page in Inkling's Biology book.)

Basic design considerations
Basic layout: The page can be read in both landscape and portrait mode. The layout ensures that the main window is the same width in each mode, so that the formatting stays the same. Portrait mode gives you reading without the clutter, landscape mode gives you access to options on the side. Should the sidebar be on the left? Should the sidebar scroll independently, and how do I do that? When publishing Our Choice, Pushpoppress decided on landscape mode only, with a novel navigation, but I not sure this would look good for a real textbook.

Continuous pages: The section is not divided up in pages. This avoids problems like figure captions spreading over 2 pages, but I think I would prefer this layout for fiction books as well: scrolling down on the iPad is a more natural gesture then side-swipe in my opinion. I don't think I'll change that.

Table of content: Navigating a science textbook is very important. The menu I have on the top left, will become too big, if I include the entire TOC. Perhaps a nested drop down menu like this? Could this be made to appear on a right-click? Alternatively, I should I simply include a link to a separate page that holds the entire TOC, like iBooks? Inkling's Biology book has a cool side-swipe feature where the TOC appears, but I have no idea how to do that.

Reading progress: Almost all readers have some indication of far along your are in the book or chapter, usually with dots. Much of this is a legacy from novels. In a textbook, you don't really care how much of the book you have read. Perhaps it's important for a chapter? Certainly in a section the position of the scroll bar is indication enough? I don't think this has high priority.

Future plans
My plan for version 2 is to add interactive figures (should they be embedded or in separate pages?), and in subsequent versions I plan to add additional items to the right side-bar. What should these items be? Movies? Quizzes? Again, suggestions are very welcome, as are links to similar mock-ups.

Now is the time to decide what interactive chemistry e-books should look like. What do you think?

Related blogposts:
Interactive chemistry ebooks: interactive figures
Interactive chemistry ebooks: highligt and annotate
iPad: even 3D molecules that can be viewed from any angle
The Molecular Modeling Basics Electronic Color Supplement

Discussion page for Section 3.1 The electron density

2011-05-05T08:25:00.000+02:00

The purpose of this post is to provide a forum for discussion of Section 3.1 The electron density of Molecular Modeling Basics. Please post your questions an comments about electron density below.

A Jmol atom picker

2011-04-04T08:41:00.000+02:00

On this page you can find a simple example of how to identify atoms in a Jmol molecule. My idea is to use it to create interactive practice problems for identifying chiral centers and things like that. But I haven't really had time to pursue it, and for a while I even forgot the link to the page. So, I'll throw it up here so I know where it is, and perhaps someone else can make use of it too.

You can see the whole code on the page (right-click, "view page source", but the main new thing here is the

jmolGetPropertyAsArray

function. The code is basically and extract of this page, written by Angel Herraez (thanks, Angel, for pointing me to that page).

Building molecules: What's in a name?

2011-03-20T10:44:00.000+01:00

In the discussion of a previous blog post Geoff Hutchison and Marcus Hanwell made me aware of a really cool building feature in Avogadro: building by naming. I demonstrate this feature on the screencast above. All the hours spend learning organic nomenclature is finally coming in useful!

Notice that, just like in MolGrabber, you get 2D coordinates, so it is important to energy minimize the structure.

ChemDoodle Web Components: 2D to 3D and MolGrabber

2011-03-13T09:41:00.003+01:00

Kevin Theisen and his colleagues over at iChemLabs have made a very useful web page that is both a great builder and a great educational tool. I have written about "2D to 3D building" before, but what makes this site special is the integration with MolGrabber combined with Chemical Indentifier Resolver by Markus Sitzmann, which generates 3D coordinates.

First, I can't think of an easier way* to build (i.e. generate the coordinates of), say, aspirin, than typing "aspirin". And of course you can modify the structure further, using aspirin as a starting point. In the screencast I show how to save the file and load it into Avogadro for minimization and GAMESS input file generation. UPDATE: Kevin has informed me that the coordinates you get are the 2D coordinates from the sketcher. *Also, there is actually an easier way as explained in the comments.

Second, the site is a great tool for showing/learning the connection been nomenclature and structure (what's the difference between 1-butene and 2-butene?), as well as the connection between 2D and 3D structure (cyclohexane is not flat like benzene!).

Currently, in order to see the 3D model you need to use the Chrome or Firefox 4 browser. Also, it is a lot of fun to draw the 2D molecules on the iPad! But the 3D model does not work in mobile Safari yet.

Lastly, iChemLabs has also made another very useful site where you can go directly from the name to the 3D structure (though here you don't have access to the coordinates).

The Molecular Modeling Basics Electronic Color Supplement

2011-03-05T12:09:00.003+01:00

Figure 3.5. (a) RHF/6-31G(d) 0.002 au isodensity surface with superimposed electrostatic potential for (a) cis-HO(H)C=C(H)OH and (b) cis-CH₃(H)C=C(H)CH₃ and. In both cases, the maximum potential value is 0.05 au.

From page 79 of Molecular Modeling Basics: "Figure 3.5 shows such plots for cis-CH₃(H)C=C(H)CH₃ and cis-HO(H)C=C(H)OH and clearly shows the difference in polarity between hydroxyl and methyl groups."

Not really.

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

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

I have now made an electronic color supplement in epub format (dowload here), which is an edited compilation of these blog posts. I see it as the next step in the evolution of Molecular Modeling Basics, and it is my first experiment with the epub format. I used the (open source!) Sigil software, as suggested by Henry Rzepa whose ebooks and wiki how-to served as an inspiration.

ePub standards are far from universally accepted, so this color supplement can look very different in different readers. On my Mac it looks OK in the Firefox epub reader plugin but not in Adobe Digital Editions or Stanza, while on my iPad it looks OK on both Stanza and iBook. I would be interested to hear about problems or successes with these and similar platforms and software.

Almost all figures link to web pages with interactive Jmol versions of the figures. None of the interactive models will work on the iPad, due to its lacks of Java support. In the not too distant future, Mobile Safari will support WebGL, and my plan is to slowly convert the interactive figures from Jmol to ChemDoodle Web Components, as the capabilities of the latter software increases.

Right now, accessing the interactive figures and videos switches you from the reader to a browser.
When the new epub format, epub3, matures (along with the readers) I hope it will become possible to view and interact with these interactive features directly within the reader.

In fact I hope to use this color supplement as a "laboratory" to experiment with these new capabilities as they become available, and make this color supplement a prototype for the next generation of scientific book publishing. Wish me luck ...

Related blog posts:
ChemDoodling: now in 3D, but not (yet) on the iPad
ChemDoodling on the iPad and the future of interactive chemistry text books
iPad: even 3D molecules that can be viewed from any angle

ChemDoodling: now in 3D, but not (yet) on the iPad

2011-02-06T10:07:00.003+01:00

You need Google Chrome or Firefox 4 to view the molecule.

In a previous post I wrote "The 3D version of ChemDoodle Web Components requires something called WebGL, which is not available in standard browsers yet, but should be soon. You can get access to it now by downloading Google Chrome (BETA)."

Now you just need Google Chrome or Firefox, i.e. if you view this post in the Google Chrome or Firefox 4 browser you should be able to interact with the molecule. Other browsers should include WebGL soon, and when it reaches Mobile Safari, it should work on the iPad.

This means ChemDoodle Web Components is ready for use in research and education. It's quite easy to use: see here for installation instructions and here for a simple html example. Furthermore, "Over the next couple years, we intend to match all functionality in JMol" according to Kevin Theisen.

I note that the ChemDoodle Web Components are open source, so feel free to pitch in!

Course: QM/MM modeling of enzymatic reactions

2011-02-06T09:23:00.000+01:00

source

I am teaching a graduate course QM/MM modeling of enzymatic reactions. Though this is about as non-basic as you can get in molecular modeling it may be of interest to some readers of this blog.

In addition to the course website, which includes a list of papers and lectures slides updated weekly, I also hope to make blog posts summarizing the in-class discussion. This will happen over at Proteins and Wave Functions, a new group blog I have started to facilitate group communication - both within the group and with the scientific community at large (constructive comments from anyone are greatly appreciated).

Jean-Claude Bradley's links to e-resources for organic chemistry

2011-01-30T11:14:00.001+01:00

An interesting e-resource page for undergraduate organic chemistry compiled by open science pioneer Jean-Claude Bradley.

I haven't gone through all the links yet, but my favorite so far is the chirality relationships exercise (number 10 on the list), where Jmol is used to ask chirality-related questions (NB: I could get this to load properly in Safari, but not Firefox or Chrome).

I think this mirrors the "mental rotations", that experienced organic chemists do to answer such questions, quite accurately. I'm not an experienced organic chemist, so I hope I didn't make too many mistakes in the screencast.

ChemDoodling on the iPad and the future of interactive chemistry textbooks

2010-12-19T13:15:00.005+01:00

Interactive model of morphine

I have bought an iPad! In honor of this purchase I bring to you this blog's first interactive figure that also works on the iPad (and most other mobile devices). It is made with ChemDoodle Web Components (a modified version of this page), an open source javascript based toolkit for chemistry, made by Kevin Theisen and co-workers at his company iChemLabs.

Readers of this blog will know that I am quite fond of Jmol for interactive molecular models, but Jmol is written in Java, which is not supported by the iOS operating system that iPads, iPhones, and iPods use - and perhaps it never will be. This decision by Apple basically means back to square one for interactive chemistry when it comes to the iPad.

I know of just two options for interactive models for the iPad: the Molecules app by Brad Larson and ChemDoodle Web Components (TwirlyMol does not appear to be interactive on the iPad). The Molecules app looks a bit more three dimensional, but works only on the iPad. Chemdoodle Web Components should work on most browsers and most operating systems, and a fully 3D version is also available. The 3D version of ChemDoodle Web Components requires something called WebGL, which is not available in standard browsers yet, but should be soon. You can get access to it now by downloading Google Chrome (BETA).

It is the Molecules app that is used in the interactive text book from Inkling that I wrote about earlier (thanks again to Henry Rzepa for the info). But I think ChemDoodle Web Components holds tremendous promise for interactive chemistry textbooks when combined with another new innovation on the horizon: EPUB3.

Epub is basically code that makes XHTML look nice when viewed in an epub reader (such as iBooks), but the current version does not allow for things like javascript, needed for interactivity. That will change with epub3 and, when combined with ChemDoodle Web Components, should allow us to make interactive chemistry textbooks that can be read on most devices. It will be an exciting time.

Computational chemistry exercises

2010-12-11T11:20:00.008+01:00

Someone (I'll call him N. O'Boyle ... no, too obvious ... Noel O.) wrote me asking if I had any computational exercises I'd be willing to share. Since I took the trouble of writing him back, it occurred to me that I had free blog material. Here's my reply in a slightly edited form.

When I taught a computational chemistry course in Iowa I used exercises from "A Laboratory Book of Computational Organic Chemistry" by Warren Hehre et al., and Spartan in thosepre-Avogadro days.Specifically experiments 4, 11, 34 and 76. See here for an example. The computational component of the other half of the course was individual research projects.

I wouldassign an experiment on a Monday, discuss it the following Monday, and have awrite-up due the Monday after that. I graded the first write-up (of exp 4) very lightly and then gave the student this example write-up of exp 4 so they could see how how to do it. I also made this check list for a report and a list of questions for each experiment (taken from the book): exp 4, exp 11, exp 34, and exp 76.

Here in Copenhagen I co-teach a similar course with 5 other people, so I just get1-2 exercises a year, and here I try to fit the content of the exercises inwith the other instructors and the topic I cover. I teach the chapter onDFT and here I have developed an exercise using bond energies.

Sometimes I also teach the chapter on geometry optimization andthen I use exp 76. Other instructors use Gaussian/Gaussviewso that's what the student tend to use here too, so I have made notutorials to go with the exercises.

"Exercises" in Molecular Modeling Basics: In Chapter 4 I illustrate applications of QM to variouschemical problems, and Chapter 5 gives you some details of the underlyingGAMESS input and output files (and there are now several blog posts witheven more information on the various examples). The intent is that peoplecan reproduce the results I present in Chapter 4 relatively easily. Depending on the level of the course, reproducing these example may be challenging enough. Otherwise, onecould easily come up with additional related problems. Let me know if thatis of interest.

Many of the molecules and concepts are very P-chem oriented, i.e. usessmall non-organic molecules to illustrate P-chem concepts, but there aresome organic molecule/concept examples too: steric strain, hydrogen bonding, amide hydrolysis.

Arsenic and odd life

2010-12-05T09:13:00.004+01:00

A recent Science paper describes "A bacterium that can grow by using arsenic instead of phosphorus", and molecular visualization is used to elegantly illustrate the basic idea, a shown in the above screencast (for the whole video visit gizmodo).

Whether the basic idea is correct is a whole other matter: for examples see here and here (be sure to check out the comments).

Simulations in teaching physical chemistry: thermodynamics and statistical mechanics

2010-12-04T12:57:00.007+01:00

In this post I summarize the simulations and I have used in teaching thermo and stat mech, and talk a bit about how I use them.

I co-teach two quite similar courses on this topic: one for nano-students and another for chemistry and biochemistry students. In the nano course we use the book Molecular Driving Forces by Dill and Bromberg, and in the other Quanta, Matter, and Change by Atkins, de Paula, and Friedman. At the end of this post I have organized the simulations by chapter for each book.

Some of simulations I have made (or modified extensively) and most of these have been discussed in previous blog posts, so I simply give the link to the respective blog post where there is more information.

The other simulations are from the Molecular Workbench (MW) library of models, and here I provide links that will open in MW, so you need to install MW before clicking on the links. For some of them I also provide a brief description of what concepts try to demonstrate using the simulations.

How do I use the simulations?
All simulations are used during lecture to visualize concepts, start discussions, and motivate equations. I'll take Illustrating energy states as an example: instead of saying "Molecules in a gas translate, rotate, vibrate, and ....", I say "Here is a zoomed-in view of butane gas where you can see the molecules. You can see that individual molecules move differently. How do they move differently? Anyone? Right, they have different speeds. This kind of motion is called translation. What else? ..."

Practical tips
On a very practical note, my own simulations are all on web sites and I make sure to open all of them before the lecture, while I have all the MW simulations for the course indexed on a single MW page (click here to open in MW). It is not possible to embed these simulations in Powerpoint slides, but you can switch between Powerpoint and other applications without quitting Powerpoint (on a Mac you use command-tab and on Windows i believe it is windowskey-tab). Note that you need access to the internet in the lecture room.

While I have screencasts of most of simulations on the blog posts, I don't use these during lecture. I think it is too passive, and puts the students to sleep. But I believe the screencasts are a good way for the students to review the main points of simulations after the lecture. I put links to the blog posts on the course web site and in the lecture notes.

Is using simulations a good idea?
If possible I try to use a simulation within the first five minutes of a lecture, and have a maximum of 20 minutes between simulations. I now only have one (45 minute) lecture left where I don't use a single simulation and I can just feel how I loose the student's attention after about 30 minutes. You can just see it. That being said, no one has ever mentioned the simulations in their course evaluations (good or bad), so I have no hard evidence that it improves my teaching. But I can tell you that I enjoy lecturing much more with the simulations, so unless I get complaints I'll keep doing it.

Making room for simulations in the lecture
I have taught the topics for many years without any simulations, and was never at a loss for material to cover. Lecture time is precious, and these simulations take time to present and discuss. You really have to introduce the simulation carefully (don't rush this part!) before you start them, and very often you want the students to speculate about what will happen before you start them. Furthermore, they tend to stimulate many more questions, that you can hopefully turn into a discussion instead of simply answering them, than derivations - that's the whole point.

So how do you "make room" for the simulations? I have cut out most of the derivations from the lectures. To pay for my sins, I provide relatively detailed (typed) lecture notes ahead of lecture (I generally don't use Powerpoint), which include step-by-step derivations. So I'll say things like "Starting with these assumptions we can write down this equation. This can be rewritten as this equation, which is much simpler. The details on how we got from here to there are in your notes, but note that in step 3 we assume that ... which is an approximation." No complaints so far. If only more progress had been made on simulating derivations ...

Here are the simulations organized by chapter

Molecular Driving Forces by Dill and Bromberg (1st edition)

Ch 6: Entropy and the Boltzmann distribution law
Illustrating entropy

Ch 10: Boltzmann distribution law
Polymer unfolding: The book uses two simple bead models of polymers in this chapter to illustrate micro and macrostates and model protein melting. I use this example extensively both in lectures and homework problems. So I made this simulation to illustrate how higher energy macrostates become more likely at higher temperatures.

Ch 11: Statistical mechanics of simple gasses and solids
Illustrating energy states
Energy states in the water molecule: a slightly more complicated molecule than HCl (used in Illustrating energy states) with more than one vibrational mode and 3 rotational degrees of freedom.
Internal energy and molecular motion
Entropy, volume, and temperature

Ch 12: Temperature, heat capacity
The molecular basis of differential scanning calorimetry: heat capacity and energy fluctuations

Ch 13: Chemical equilibria
Seeing chemical equilibrium (opens in MW)
Dalton's law of partial pressure (opens in MW)

Ch 14: Equilibria between solids, liquids, and gasses
Seeing specific and latent heat (opens in MW): I use this simulation to illustrate how the same substance can be solid, liquid, and gas depending on the temperature.
A gas under a piston (opens in MW): I use this simulation to show that, for example, decreasing the pressure can have the same effect as increasing the temperature.
The phase diagram explorer (opens in MW)
Raoult's law: ideal solutions (opens in MW): Here, I use the simulation of the pure liquid to illustrate vapor pressure.

Ch 15: Solution and Mixtures
Mixing gasses, and mixing of ideal and non-ideal liquids
Raoult's law: ideal solutions (opens in MW)
Raoult's law: negative deviation (opens in MW)
Raoult's law: positive deviation (opens in MW)

Ch 16: Solvation and transfers of molecules between phases
Visualizing osmotic pressure in an osmotic equilibrium (opens in MW)
Desalination using reverse osmosis (opens in MW)

Quanta, Matter, and Change by Atkins, de Paula and Friedman (1st edition)

Ch 13: The Boltzmann distribution
Illustrating energy states
Energy states in the water molecule: a slightly more complicated molecule than HCl (used in Illustrating energy states) with more than one vibrational mode and 3 rotational degrees of freedom.
Internal energy and molecular motion

Ch 14: The first law of thermodynamics#
The molecular basis of differential scanning calorimetry: heat capacity and energy fluctuations

Ch 15: The second law of thermodynamics
Illustrating entropy
Entropy, volume, and temperature

Ch 16: Physical equilibria
Seeing specific and latent heat (opens in MW): I use this simulation to illustrate how the same substance can be solid, liquid, and gas depending on the temperature.

A gas under a piston (opens in MW): I use this simulation to show that, for example, decreasing the pressure can have the same effect as increasing the temperature.

The phase diagram explorer (opens in MW)
Raoult's law: ideal solutions (opens in MW): Here, I use the simulation of the pure liquid to illustrate vapor pressure.

Raoult's law: negative deviation (opens in MW)
Raoult's law: positive deviation (opens in MW)

Mixing gasses, and mixing of ideal and non-ideal liquids

Visualizing osmotic pressure in an osmotic equilibrium (opens in MW)
Desalination using reverse osmosis (opens in MW)

Ch 17: Chemical equilibria#
Seeing chemical equilibrium (opens in MW)
Dalton's law of partial pressure (opens in MW)

# I don't teach this part of the course, but if I did I would use these simulations

Related posts:
An Atkins Diet of Molecular Workbench
One, Two, Three, MD
Tunneling and STM (a first stab at using Molecular Workbench to teach quantum mechanics)

Illustrating mixing

2010-12-04T12:48:00.008+01:00

This screencast shows Molecular Workbench simulations I have made to illustrate mixing.

The first simulation illustrates the mixing of 2 ideal gases, which mix readily. Since the gas particles don't interact you can think of the mixing as each gas expanding to fill both containers independently of each other. As I have shown in this simulation, the driving force for this expansion is an increase in entropy. Therefore, the driving force for mixing two ideal gasses is also purely entropic.

The second set of simulations illustrates the mixing of 2 liquids. Since they are liquids there must be attractive interactions between the atoms. If there were no interactions they would be (ideal) gasses. The strength of the interactions (and the temperature) determine whether they mix or not.

In the first liquid simulation, the attraction between two green atoms (ε_GG), between two blue atoms (ε_BB), and between a green and a blue atom (ε_GB) are the same.

This means that a green atom doesn't care whether it is sitting next to a blue atom or another green atom. The net effect is that green and red atoms are equally likely to be on the right or left side of the container, and the liquids mix for the same reason as the ideal gasses mix: the driving force is purely entropic. That means the enthalpy of mixing is zero:

This is the definition of an ideal mixture (or ideal solution). The two liquids will mix at any temperature.

In the second liquid simulation, the attraction between two blue atoms (ε_BB) is stronger than between two green atoms (ε_GG) and between a green and a blue atom (ε_GB).

Note that the ε's are negative: a smaller ε means a stronger attraction. This means that the blue particles would rather be with other blue particles, i.e. the enthalpy increases if the particles are mixed.

(z is the number of contacts between particles in solution, and x_G is the mole fraction of green atoms). This is an example of a non-ideal mixture, where the definition for a non-ideal mixture is

Because Δ_mixH > 0 this non-ideal solution mixes spontaneously (i.e. Δ_mixG < 0) only for

Oil and water is a common example of such an non-ideal mixture: the oil-oil interactions are stronger than the oil-water and water-water interactions.

Salt and water is another example of on idea mixture, but here Δ_mixH < 0 so salt and water almost always mixes spontaneously. The interpretation is that the interactions between the salt ions and water is stronger than the average interaction between salt ions and between water molecules.

Implications and limitations
The definition of Δ_mixH in terms of the ε's suggest that liquids should also mix if

which would be a more general definition of an ideal mixture.

This is tested in the third liquid simulation. As you can see the liquids mix more than in the second simulation, but not quite as much as in the first simulation. This is mostly because of the simulation runs only for 100 picoseconds, which is to short to mix fully. But another reason is that there is less space (on average) between the blue atoms compared to the green atoms, because the blue atoms attract each other more. The next effect is that the blue particles tend to stay together to lower the enthalpy. More mathematically, z (the number of contacts between particles in solution) is not exactly the same for the blue and green particles so the interpretation of Δ_mixH in terms of the ε's breaks down. The "safest" definition of an ideal mixture thus remains:

i.e. "like dissolves like".

Accessing the simulations
You can play around with the simulations here and here, or you can download the models here and here if you have Molecular Workbench installed on your computer.

The Open Science unjournal

2010-11-21T16:15:00.014+01:00

In this post I do something a little different than normal for this blog: I introduce a scientific unjournal called Open Science. It doesn't exist, but it should.

What is an unjournal?
An unjournal is to journals what an unconference is to conferences. To define what an unjournal is, take the first 2 sentences in the wikipedia entry on unconferences and substitute a few words: "An unjournal is a facilitated, participant-driven journal centered on a theme or purpose. The term "unjournal" has been applied, or self-applied, to a wide range of publications that try to avoid one or more aspects of a conventional publishing, such as loss of copyright, high fees, [list your favorite pet-publishing-peeve here]."

Here are some frequently asked questions (FAQs) about Open Science:

How do I publish in Open Science?
The usual way is to deposit your manuscript on openscienceunjournal.org. Once you have completed the submission process the paper is given a time and date stamp and the paper is published and open for review (see next question). While there is a recommended template available for the paper, there is no fixed format. It is possible to upload all your raw data or link to the data if it is hosted elsewhere. Experience has shown that this tends to increase the scores (see below) of your papers significantly.

You can also choose to publish the paper on any of the unjournal sites whose contents are linked to openscienceunjournal.org. These sites are often run by established publishers and offer more user-friendly interfaces, but may require a fee and may ask you to give up your copyright (though the content is open access by definition).

Is Open Science peer reviewed?
For a publication in an unjournal such as Open Science the question should be: is my article in Open Science peer reviewed? That is in large part up to you. The paper is open for online, non-anonymous, and completely transparent review and you have 2 months in which you can change the content of the article in response to the comments. After the 2 month period you cannot change the content, but you can of course respond to new comments online. For very serious criticisms you may want published a new Open Science article to respond.

It is up to you to solicit reviews, though any published author of a peer-reviewed paper (defined below) can review your paper during the 2 months review process. Based on our experience only well-written and well-presented articles on scientifically interesting questions get reviewed.

An Open Science article is neither rejected nor accepted at the end of the review process. Instead it receives an initial score (see next FAQ). Obviously, any paper that does not generate a single review (or is reviewed but gets an initial score of 0) is not considered peer reviewed.

What is the impact factor of Open Science?
For a publication in an unjournal such as Open Science the question should be: what is the impact factor of my article in Open Science? During the 2 months review process each reviewer gives the paper a score between 5 (good) and 0 (bad). This score can be adjusted by the reviewers based on the changes you make within the first 2 months, after that it is fixed. The initial score for your paper is the average of all reviewers final scores.

If the paper is cited it receives an additional score (called the current score). The score is determined by the number of citations, and if the citing paper is published in Open Science, the score is weighted by the initial and current score of that paper. The authors of that paper can also choose to indicate how important your paper was to theirs. If high scoring papers cite your paper in a positive manner, the current score of your paper increases. Self-citations are not included.

Why should I review for Open Science?
The work you put into reviewing is now documented for all to see. Have you contributed greatly to science by identifying Open Science papers with high current scores? Do the reviews you submit carry more weight with the author and other reviewers as a result? Some sites now list Open Science reviewers with particularly high impact as a kind of editorial board for the journal.

How should I cite an Open Science paper?
One suggestion is: Author(s), Title, Open Science, date of submission, initial/current score. If you publish in Open Science using the suggested template, the current score is updated automatically.

Why should I publish in Open Science?
There are many reasons:
(1) You retain the copyright and anyone can see the paper.
(2) Your paper is accessible upon submission. (Don't rush to publish though: you only have 2 months to get a good initial score).
(3) The impact of your paper is evident in the citation, but disconnected from the conventional impact factor of the journal you managed to get it in to. The initial score of your paper can help the paper off to a good start, but your truly important papers will ultimately be identified by its current score.
(4) You choose the publishing format you like. What's your pleasure? machine readable? interactive figures? link to raw data?
(5) Your paper is a living document: comments or questions continue to roll in on important papers and you can update links to your papers (related articles, a new data format) as you see fit.
(6) If you write a good paper, you will get more reviews (i.e. more suggestions and input) but the rantings of a single idiot reviewer will not prevent publication. Isn't this the place to publish daring and ground-breaking work?

So what brought this on?
The blogosphere: Egon Willighagen's latest post got me thinking about this particular idea, but the general problems it is trying to address was brought to my attention by many other blogs such as Michael Nielsen's The Future of Science and Is Scientific Publishing About to be Disrupted? posts; most posts by Peter Murray-Rust; Henry Rzepa's long fight to include interactive figures in conventional journals; Mat Todd's excitement for an unconference and the discussion it generated at Derek Lowe's blog. Why can't we have this in a journal?

Is the Open Science unjournal a good idea?
The blogosphere will decide: no comments on, and no re-blogging of, this post will mean this idea dies a quiet death (by receiving an initial and current score of 0). But if you get enough smart people fired up about an important idea that can be solved by IT, good things can happen.

The molecular basis of differential scanning calorimetry: heat capacity and energy fluctuations

2010-11-06T21:14:00.003+01:00

Melting and boiling points are convenient and important measures of stability. But how do you measure a melting point of, for example, nanoparticles that are too small to see?

This screencast shows two Molecular Workbench simulations (you can find them here and here) I made [see credits at the end of the post] to illustrate the connection between phase transitions, changes in heat capacity, and energy fluctuations, and the slides below takes you through the basic ideas behind them.

Slide 1: In the first simulation heat is added to a nano-particle and the resulting temperature increase is measured. When viewing the simulation notice that the temperature increases less during the melting/evaporation.

Slide 2: To analyze the data we first switch the x and y-axes, so that heat added (i.e. the internal energy, U) is plotted as a function of temperature.

Slide 3: The data is a bit noisy (mainly because the simulation heats the particle too fast: from 0 to almost 3000 K in about 200 picoseconds!), so I smooth it by fitting a curve to it.

Slide 4: From the smoothed data I can calculate how fast the energy changes with temperature. This is the heat capacity (C_v), which peaks at a temperature around 1350 K - the melting temperature of the particle.

Slide 5: This observation forms the basis of differential scanning calorimetry, which measures the temperature as a function of the flow of energy to a system, and determines the melting point by finding the temperature where the heat capacity peaks.

Slide 6: One way to explain why the heat capacity peaks at a phase transition such as melting is through its relation to energy fluctuations: the system changes most during a phase transition ("bonds" between particles are broken and formed), so the energy fluctuates more, meaning that the heat capacity is largest.

Slides 7, 8, and 9: In the second set of simulations the energy is plotted a function of time at 3 temperatures: before the particle melts (500 K), when the particle melts/evaporates (1350 K), and after the particle has evaporated (3000 K).

Slide 10: Results from the 3 simulations are compared. Clearly the fluctuations are largest when the temperature is 1350 K. The fluctuations at 3000 K are larger than at 500 K, even though the heat capacities are similar. This is because the heat capacity is proportional to the average energy fluctuation divided by the temperature squared (slide 6).

Calorimetry

View more presentations from molmodbasics.

You may wonder why we don't see two heat capacity peaks: one for melting and one for evaporation. This is because of the particle is so small (i.e. composed of relative few particles). For a macroscopic systems (like water) the phase transitions are well defined. Water is ice at 272 K, melts at 273 K, and is a liquid at 274 K (at 1 atmosphere of pressure); and the heat capacity has a very narrow peak at 273 K.

As particles become smaller their phase transitions become less well defined, the heat capacity peak becomes broad, and in some cases (like this one) you get a single heat capacity peak for melting and evaporation. This means that the phase transition cannot really be classified as melting or evaporation and that is occurs over a relatively large temperature range. Li and Truhlar have an interesting article on this subject.

If you would like to play around with or modify the simulations they can be downloaded here and here, but you need to download Molecular Workbench first.

Credits: The simulations are based on models and scripts by Arie Aziman and Carlos Gardena, who based their work on a model by Dan Damelin, i.e. they are made possible, like Molecular Workbench itself, by open source science.

Nobel reactions

2010-10-07T20:40:00.003+02:00

You may have heard this already: this years Nobel prize in chemistry went to Richard F. Heck, Ei-ichi Negishi, and Akira Suzuki "for palladium-catalyzed cross couplings in organic synthesis". Like last year, molecular modeling and animation is here to bring the science behind the prize to life. This time, through the very excellent Chemtube3D (a Molecular Modeling Basics favorite) and its page on on organopalladium chemistry.

Illustrating energy states

2010-10-02T13:38:00.001+02:00

Here is a screencast showing 2 simulations I made to illustrate energy states (often also called microstates). Open any book on statistical mechanics and you'll see a formula like this (for the Boltzmann distribution),

and, perhaps, a figure much like this.

But what are these "energy states" and ε's exactly? Hopefully the screencast and associated simulations, which you can find here and here, will help make this clearer.

About the simulations
If you want to play around with the HCl simulation, note that you have to reload the page if you want to change the translational motion. This simulation is made with Jmol. The molecular dynamics simulation of butane is made with Molecular Workbench.

iPad: even 3-D molecules that can be viewed from any angle

2010-09-11T12:43:00.003+02:00

I recently came across this report on inkling.com, a new company aimed at bringing interactive textbooks to the iPad. Interactive molecular models was mentioned two times and clearly left an impression on the author (italics are mine):

"Inkling’s software turns textbooks into interactive content, with video, hyperlinks between text and images, notes that can be shared between students and teachers, and even 3-D molecules that can be viewed from any angle."

"MacInnis – who worked at Apple for eight years, including a stint in the company’s educational division — says that the iPad is the perfect device for the kind of interactivity that Inkling provides because it has the ability to produce high-end graphics, such as the 3-D spinning molecule that is a feature of the company’s biology textbook."

This feature is also shown in inklings promotion video, excerpt below:
Update: MacInChemBlog keeps a list of science related iPhone/iPad apps.

Molecular Modeling Basics reviewed in Chemistry World

2010-09-02T20:11:00.001+02:00

A very nice review of the book (and blog!) by fellow blogger Henry Rzepa in Chemistry World. Much appreciated.

Entropy, volume and temperature

2010-08-29T19:11:00.005+02:00

This screencast shows a Molecular Workbench simulation I made to illustrate the connection between entropy, volume, and temperature.

Each simulation runs for 50 picoseconds (ps) and records how much time the two particles spend together (as a dimer) and apart (as monomers). These times are a reflection of the relative probabilities of the dimer (p_dimer) and monomers (p_monomers).

In the first simulation of the screencast the particles spend 38.3 ps together and 11.7 ps apart. When I double the volume (while keeping the temperature constant) in the second simulation, the particles spend 32.7 ps together and 17.3 ps apart.

This makes intuitive sense: when the particles are not together they have a harder time finding each other again in a bigger volume. Put another way, the probability is lower that the particles find each other when they move in a larger volume.

Now I want to connect this intuitive explanation with a thermodynamic one (i.e. an explanation in terms of free energy and entropy) :

The relative probability of being together and apart is given by the change in free energy (ΔA) on going from the dimer to the two monomers (dimer -> 2 monomers)

The increase in p_monomers/p_dimer (from 0.31 to 0.53) on doubling the volume must mean that ΔA decreases when the volume is doubled. Statistical mechanics tells us that the only thermodynamic term in A that is affected by volume (V) is the translational entropy

and that an increase in volume will increase the entropy of a particle

So when the volume is doubled the entropy of each particle (the dimer and each monomer) is increased by Rln(2). As a result ΔS for the reaction dimer -> 2 monomers increases [by Rln(2)] and ΔA = ΔU - TΔS decreases by RTln(2). So the probability of the dimer relative to the monomers decrease because the entropy is increased when the volume increases.

In the last simulation of the screencast I double the temperature (while keeping the volume constant) compared to the first simulation. As a result the particles now spend less time together (19.2 ps) than apart (30.8 ps).

This also makes intuitive sense: the dimer is more likely to be struck by a particle with a kinetic energy larger than the potential energy that is holding it together: Furthermore, when the monomers happen to collide they are more likely to have a combined kinetic that is larger than the potential energy holding the dimer together, i.e. "with too much kinetic energy to stick together".

As before, this must mean that ΔA decreases when the temperature is increased. The increase in translational entropy due to temperature is indeed larger than for the volume

and consistent with the observed larger change in the time the particles spend apart.

If you think of the dimer as a crude model of a salt you want to dissolve, this explains why dilution (increasing the volume) or heating increases the solubility (the amount you can dissolve). Conversely, increasing the concentration (the amount of particles per volume) or decreasing the temperature helps dimer formation and, in a larger sense, self assembly.

You can play around with it on this web page, or you can download the model if you have Molecular Workbench installed on your computer. Enjoy!

Clarifications and approximations
You may wonder why I discuss the Helmholtz free energy change (ΔA; some books use ΔF) rather than the Gibbs free energy change, when the volume changes. This is because the volume of the system (i.e. the two compartments) does not change. Another argument is that no work is done during the expansion.

I have implicitly invoked the ergodic hypothesis which states that the probability computed using a collection of molecules at a single instant in time is equal to the probability computed for a single molecule over a long period of time. While it makes intuitive sense, I don't believe this has been rigorously proven.

Related blog posts
Illustrating entropy
Where does the ln come from in S = k ln(W)?

Installing GAMESS on a linux PC

2010-08-24T08:38:00.007+02:00

Probably best viewed in full screen mode

Alchemist requested a post on installing GAMESS on a linux system: after all "I[f] someone can not install a program, how he/she can use it?" Fair enough.

If you have a mac or windows machine, pre-compiled versions of GAMESS is available (see this post on for installing such a version on Macs). But if you have bought a linux computer, you do have to compile GAMESS yourself.

Three disclaimers:
(1) The post is aimed at people who want to install GAMESS on their personal computer that happens to run linux.
(2) The post is based on installing the March 25, 2010 version of GAMESS and corresponding scripts as they were distributed in mid August, 2010. If you view this post a few years from now things may have changed.
(3) I don't have a linux PC, so the screencast is made on a shared research computer. This affects one of the steps as I describe below, and may affect others steps that I don't know about.

The screencast assumes you have already downloaded gamess (gamess-current.tar.Z). To download GAMESS start here and follow the instructions. If you have problems watch the first few minutes of the screencast this post.

Here's a summary of the steps in the screencast:

1. Unpack GAMESS: type "zcat gamess-current.tar.Z | tar -xvf -"
Now you should have a GAMESS directory. The file gamess/misc/readme.unix has most of the instructions you need. What follows it basically a condensed version with a few modifications.

2. In the GAMES directory: type "./config".
Follow the instructions and answer the questions. You need to know whether you have a 32- or 64-bit chip.

3. Type "/sbin/sysctl -w kernel.shmmax=1610612736"
If you are curious about what this does, read section 5 of the file gamess/ddi/readme.ddi.
You need to be logged in as root for this to work. I don't have root access on the machine I used, so I get an error message in the screencast. But our system administrator had executed the command previously.

5. Change to the ddi subdirectory (gamess/ddi). Type "./compddi >& compddi.log" and then "mv ddikick.x .."

6. Go back to the gamess directory and type "./compall >& compall.log"
It will take 30-60 minutes to compile GAMESS.

7. Type "./lked gamess 00 >& lked.log"

8. Edit the file rungms. Here are the lines I modify or add:

set SCR=./
set USERSCR=./
set GMSPATH=./

setenv ERICFMT ./ericfmt.dat
setenv MCPPATH ./mcpdata

setenv  MAKEFP $USERSCR/$JOB.efp
setenv   GAMMA $USERSCR/$JOB.gamma
setenv TRAJECT $USERSCR/$JOB.trj
setenv RESTART $USERSCR/$JOB.rst
setenv   INPUT $SCR/$JOB.F05
setenv   PUNCH $USERSCR/$JOB.dat

   if ($os == Linux)   set GMSPATH=./

#  if ($NCPUS == 1) then
#     set NNODES=1
#     set HOSTLIST=(`hostname`)
#  endif
#
   set HOSTLIST=(`hostname`:cpus=$NCPUS)
   set NNODES=1

#           echo I do not know how to run this node in parallel.
#           exit 20

9. Edit one line in the file runall: #chdir /u1/mike/gamess.
Type "./runall >& runall.log"

10. Go to the tools/checktst subdirectory and type "./checktst"
(Update: if the gamess directory is not in your home directory, you need to change this in the checktst file. See this comment)

Running GAMESS in parallel
If you computer has more than one core, you may want to run GAMESS in parallel. To run exam01 on 2 cores type "./rungms exam01 00 2 >& exam01.log &"

Other useful programs and getting started with GAMESS
Now that you have GAMESS running you may want to install other useful programs such as Avogadro, GAMESSQ, and MacMolplt.
Here are some related blogposts to wet you appetite:
Using GAMESSQ
Building a complicated molecule with Avogadro

And here are a few blogposts on getting started with GAMESS:
Some GAMESS input basics
A typical set of GAMESS calculations

Acknowledgment: Mike Schmidt and Casper Steinmann helped with this post.

August 30, 2010 update: Alchemist has made the following page on installing GAMESS on Ubuntu 64 bit. Be sure to check out the comments for additional useful links.

October 20, 2010: Are you using Ubuntu 9.10 and 10.04, Linux Mint 8 and 9, Knoppix 6.2.1 or Suse 11.2? Be sure to check out this comment by Mott.

Amide hydrolysis, revisited 2

2010-08-11T15:19:00.006+02:00

In a previous post I discussed the TSs of acetamide hydrolysis and hydrolysis of a related amide (3). I have already made a post on how to find the TS for 3, and in this post I summarize the two ways of finding the TS for acetamide hydrolysis, that I describe in Chapter 5 of the book (where you can find many more details).

To start, I use Avogadro to build the geometry shown in Figure 5.32a and use it to construct the input file for an optimization with 4 constrained bond length. The values for the constraints are taken from a paper. Figure 5.32b shows the equilibrium geometry obtained with GAMESS, after 17 steps. Based on this geometry GAMESS finds the TS in three steps.

Figure 5.32a. Initial guess geometry for the constrained optimization discussed in Figure 5.33.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

Figure 5.32b. The geometry resulting from the constrained optimization and the normal mode of the imaginary frequency computed for this structure at the PM3 level.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

It will not always be possible to find an article that reports the TS structure for your reaction of interest. So let’s try to find the TS without the information from the article. I start from the same structure (Figure 5.32a), and I arbitrarily pick C1-O10 as the reaction coordinate constrained to 1.70 Å.

This constrained optimization results (after 23 steps) in the geometry shown in Figure 5.35a. The subsequent TS search results in the geometry shown in Figure 5.35b, which has no imaginary frequency and is clearly the product.

Figure 5.35a. The geometry resulting from a constrained optimization in which the distance between atoms 1 and 10 (see Figure 5.32a for numbering) was constrained to 1.7 Å, and the normal mode of the imaginary frequency computed for this structure at the PM3 level.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

Figure 5.35b. The geometry resulting from a TS search initiated from the geometry shown in Figure 5.35(a).

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

The fact that the C–O bond is re-formed indicates that it should be stretched more in the initial guess, so I repeat the optimization with the C–O distance constrained to 2.0 Å instead of 1.70 Å. This results (after 39 steps) in the geometry shown in Figure 5.36. Using this as a starting geometry for the TS search leads to the TS in Figure 4.30a after 17 steps.

Figure 5.36. The geometry resulting from a constrained optimization in which the distance between atoms 1 and 10 (see Figure 5.32a for numbering) was constrained to 2.0 Å, and the normal mode of the imaginary frequency computed for this structure at the PM3 level.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

In the screencast below I try to reproduce the calculations that lead to these figures. Because I rebuild the structure in Figure 5.32a, I get different starting coordinates, so the energies and number of steps are different than what I discuss in the book for the 4-constraints approach.

In the case of the 1-constraint approach, I actually manage to find the TS using the 1.70 Å constraint. Thus, the structures in Figures 5.35b and 5.36 do not appear in the screencast. This just goes to show how finicky TS searchers are to starting geometries.

While editing the screencast I noticed that the PM3 energies of the two TS structures I find are very different (11 kcal/mol). This is also true for the TSs I found when writing the book (where they are different by 6.6 kcal/mol).

This appears to be a problem that PM3 has with structures where bonds are partially broken or formed. I could not reproduce this for structures with normal bond lengths such as the reactants and products. Note that in the book, I use M06/6-31G(d) single point energies to compute barriers, and that the M06/6-31G(d)//PM3 single point energies of the TSs found with the two different method are only different by 0.03 kcal/mol.

Amide hydrolysis, revisited

2010-07-29T20:41:00.003+02:00

Figure 4.29. Sketch of hydrolysis reaction of several amides together with their experimentally observed half-lives. The free energies of activation are computed via Equations (4.30) and (4.31). (Adapted from N. M. Hernandes et al 2008. Journal of Organic Chemistry 73: 6413–6416.)
From Molecular Modeling Basics CRC Press, 2010
January, 2011 update: The figure in the book is missing an N atom in structures 4 and 5

Figure 4.30a and b shows the TS geometries for the hydrolysis of 1 and 3 computed at the PM3 level of theory, together with the normal modes associated with the imaginary frequencies (2336i and 239i cm^–1, respectively).

Figure 4.30a. PM3 geometries of the TSs for hydrolysis of compound 1 (Figure 4.29), and the normal modes associated with the imaginary frequencies.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

Figure 4.30b. PM3 geometries of the TSs for hydrolysis of compound 3 (Figure 4.29), and the normal modes associated with the imaginary frequencies.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

As I show in the book, the predicted activation free energy of 1 is 4.2 kcal/mol higher than the activation free energy for the hydrolysis of 3. The source of the difference is roughly half electronic (1.9 kcal/mol) and half thermodynamic (2.3 kcal/mol). The explanation for the latter is the loss of translational entropy associated with hydrolysis of 1, but it is not obvious why the electronic activation energy should be lower.

For example, one might imagine that it would be energetically unfavorable to fold the chain of 3 into a ring due to some kind of strain when forming the transition state. This can be tested by studying the hydrolysis of 2 (Figure 4.29), which should have roughly the same amount of ring-strain associated with the reaction.

Indeed, the free energy of activation for hydrolysis of 2 is considerably higher than that for 3, and due entirely to an increased electronic activation energy. This points toward the importance of the amine group in the middle of the chain in lowering the electronic barrier for 3-hydrolysis, presumably by stabilizing the partially positive –NH₃⁺-like portion of the TS (Figure 4.31).

Figure 4.31. 0.002 au isodensity surface with superimposed molecular electrostatic potential of the TS for hydrolysis of 3. The maximum potential value is 0.05 au, and the level of theory is M06/6-31G(d). The orientation is the same as Figure 4.30.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

It might be tempting to ascribe the more open TS structure in 3 compared to 1 (Figure 4.30) to ring-strain, but the TS for “methanolysis” of 1 is equally open (Figure 4.32). This is presumably due to steric hindrance of the methanol methyl group and the carbonyl oxygen.

Figure 4.32. PM3 geometries of the TSs for methanolysis of compound 1 (Figure 4.29), and the normal modes associated with the imaginary frequencies.

Click on the picture for an interactive version.

From Molecular Modeling Basics CRC Press, 2010

I will discuss how to find the TS structure for the hydrolysis of 1 in a future post. I have already discussed how to find the TS structure for the hydrolysis of 3 here and here. This post describes how to verify that the TS connects the correct reactants and product, and this post describes the relationship between half lives and activation free energy.

Melting: a simple model

2010-07-25T14:26:00.002+02:00

Figure 4.26. PM3 optimized structure of the V-shaped water timer. Shown also is the normal mode corresponding to the lowest vibrational frequency (36 cm–1).
Click on the picture for an interactive version.
From Molecular Modeling Basics CRC Press, 2010

While the lowest energy conformation of three water molecules is the ring structure (Figure 4.17), there is another minimum (at least on the PM3 potential energy surface - Figure 4.26) that is 7.7 kcal/mol higher in (electronic) energy.

The free energy difference between the cyclic and V-shaped structure is zero at around 480 K. This can be considered a very simple model for melting (i.e., the T at which higher enthalpy conformations are most probable because of entropy). The entropic term has two basic contributions: there are more higher-energy structures (they are more disordered so there are more ways to make them: e.g. there are 3 identical V-shaped structures), which lowers the conformational free energy, and the structures are “floppier” (they have more low-frequency vibrational modes), which lowers the vibrational free energy.

Figure 4.28a. Structure of one of the water clusters found by Maeda and Ohno. The coordinates are taken from their supplementary materials. (From S. Maeda and K. Ohno, 2007. Journal of Physical Chemistry A. 111: 4526–4534.)
Click on the picture for an interactive version.
From Molecular Modeling Basics CRC Press, 2010

Of course, this is a hypothetical melting transition because the cyclic “ice” structure already sublimates at 285 K. The main reason for the high melting temperature is the small number of higher-enthalpy conformations, which increases quickly with the number of water molecules. For example, for water octamer [(H2O)8] a study by Maeda and Ohno found 164 different conformations, and only seven of these can be classified as some variant of the lowest-enthalpy cubic conformation (Figure 4.28a), while the rest are more disordered (e.g., Figure 4.28b). The study estimates that the temperature at which the cubic and more disordered structures become equally probable (i.e., the melting temperature) is around 280–320 K, which is significantly closer to the melting temperature of bulk ice of 273 K. The uncertainty in the estimate of Tmelt comes from the difficulty in estimating the effect of BSSE and anharmonic effects.

Figure 4.28b. Structure of one of the water clusters found by Maeda and Ohno. The coordinates are taken from their supplementary materials. (From S. Maeda and K. Ohno, 2007. Journal of Physical Chemistry A. 111: 4526–4534.)*
Click on the picture for an interactive version.
From Molecular Modeling Basics CRC Press, 2010

* The Jmol structure does not correspond exactly to the picture. When making the figure I forgot to write down which of the 164 structures I used for this figure. Still, the point is the same.

Fold.it - Solve puzzles for science

2010-06-16T21:50:00.004+02:00

I just came across this amazing game/puzzle/software - call it what you will - Foldit.

It's an interactive protein folder that evaluates the energy on the fly to let you know how well you are doing.

But, and here is the important bit, the energy function is the same used in the Rosetta packages - a state-of-the-art protein folder (at least that's how I read it - the details are a bit sketchy). This is not a toy, but an exciting experiment to test whether human intuition and pattern recognition can compete with Monte Carlo sampling.

The screencast below (see a bigger version here) shows 6 intro puzzles that teaches you some of the techniques for protein folding. You get the program here and make an account for yourself here (making an account using Foldit did not work for me).

When you start Foldit the first time it takes you into the intro puzzles immediately (i.e. bypassing the menu).

I don't know where to start about the many implications of this program for fundamental science, teaching, and public outreach. So I'll simply say wow!, and go on to level 2-3. Not a moment to lose, CASP9 is underway you know!

June 20th, 2010 update:

Not sure why I didn't think to check youtube immediately but there are many videos about foldit. Here are a couple:

Geometry and molecular motion

2010-06-07T21:59:00.005+02:00

Figure 4.18. The relative energy of H2 as a function of H–H distance and the energy of the lowest vibrational level (horizontal line).
From Molecular Modeling Basics CRC Press, 2010

When discussing structure and other molecular properties it is important to remember that the molecule is in constant internal motion due to vibrational motion. The bond length and angles that one so carefully computes and reports to so many decimal places actually changes considerably with time.

Consider the H₂ molecule. Using B3LYP/6-31G(d) and the harmonic oscillator approximation, the frequency for the H–H stretch vibration is 4450 cm^–1, which corresponds to 6.4 kcal/mol of kinetic energy. If we compute the energy as a function of the H–H bond length (Figure 4.18), we see that this kinetic energy is enough to compress the H–H bond length to roughly 0.64 Å and increase it to 0.88 Å: the classical turning points.

The classical turning points can be estimated using the harmonic approximation of the potential energy surface (PES)
where ν_i is the vibrational frequency in units of cm^–1, with a corresponding normal coordinate l_i (see section 1.3 of the book), and n is the vibrational quantum number. Thus, ν_i = 4450 cm^–1 corresponds to a turning point at l_i =±0.087 amu^½ × Å which can be converted to Cartesian coordinates (and hence a bond length) using the normal mode component L_ij In the case of H₂, this gives classical turning points at 0.621 Å and 0.865 Å, which are in good agreement with the previous estimate from Figure 4.18. The screencast below shows ho MacMolPlt can be used to estimate the classical turning points.

For lower frequencies the displacements at the classical turning points can be quite substantial, both because the PES is flatter and because higher vibrational levels are populated. For example, in the case of ethane the lowest frequency is 310 cm^–1 [at the B3LYP/6-31G(d) level of theory] and corresponds mainly to rotation about the CC bond (Figure 4.20a). For the ground vibrational level the classical turning points occur at l_i =±0.330 amu^½ × Å , which corresponds to a dihedral angle change of ±14° (Figure 4.20b, see screencast below). At room T 4% of the ethane molecules have three quanta of energy in this mode (n = 2), corresponding to a kinetic energy of 2.2 kcal/mol, enough to change the dihedral angle by ±31°.

Figure 4.20. (a) Normal mode associated with the lowest frequency computed for ethane using B3LYP/6-31G(d). (b) Structure obtained by displacing along the mode by 0.330.
Click on the picture for an interactive version.
Click here for a pop-up window
From Molecular Modeling Basics CRC Press, 2010

Some of the lowest frequencies you’ll observe are for intermolecular interactions such as the water dimer hydrogen bond (Figure 4.22). Curious how much the geometry is displaced? Why not try it yourself? (See this post to get started with GAMESS).

Figure 4.22. (a) The PM3 normal mode corresponding to (a) torsional motion about the hydrogen bond and (b) hydrogen bond stretch in the water dimer. The corresponding frequencies are 128 and 415 cm^–1.
Click on the picture for an interactive version.
Click here for a pop-up window
From Molecular Modeling Basics CRC Press, 2010

Using MacMolPlt to compute the displaced geometry
As you can see from the second equation, the normal coordinate is basically a scale factor. Thus, MacMolPlt can be used to compute a displaced geometry by converting l_i to a percent and changing the units to amu^½ × Bohr by dividing by 0.52918 Å/Bohr. Thus, l_i =±0.330 becomes 62.0% (note that there is a mistake in the book: I forgot about the unit conversion in Figure 5.27). Unfortunately, MacMolPlt does not include the m^-½ so this is only an estimate.

The Computational Microscope

2010-05-29T08:54:00.001+02:00

A great title and an interesting talk about one of the fundamental roles of molecular modeling.

Many-body effects in the water trimer

2010-05-08T13:30:00.001+02:00

Figure 4.17. M06/6-31G(d) optimized geometry of the water trimer.
Click on the picture for an interactive version
From Molecular Modeling Basics CRC Press, 2010

In a previous post I showed how molecular electrostatic potentials can be used to illustrate polarization due intermolecular interactions. Polarization has an interesting and important consequence called the many-body effect, which occurs in systems with many relatively strong interactions such as the water trimer (Figure 4.17).

The many body effect is due to the fact that polarization of the charge density of, say, water molecule A due to the A–B hydrogen bond increases the strength of the A–C hydrogen bond, which lowers the energy further. As a result the binding energy of the water trimer (i.e. its energy relative to that of three free water molecules) is lower than the sum of the three hydrogen bond strengths.

The A–B hydrogen bond strength can be gotten by (1) deleting water C and computing the energy of the A-B dimer (don't re-optimize), and (2) comparing this energy to that of two free water molecules. And similarly for the A–C and B–C dimers. The coordinates can be extracted from the Jmol interactive figure as shown in a previous post, and this post shows how to perform the calculations with GAMESS.

The increased hydrogen bond strengths are also reflected in the hydrogen bond lengths, which are 0.09Å shorter in the trimer compared to the dimer. You can measure the distances between atoms in the Jmol interactive figure by double clicking on the atoms, and water dimer geometry can be found in a previous post.

Putting it out there

2010-05-04T18:46:00.001+02:00

An inspiring post over at Noel O'Blog on sharing your PowerPoint presentations using Slideshare. Three of Noel's presentations (Quantum Chemistry I, II, and III)deserve special mention here because (a) they relate to teaching molecular modeling and (b) the third mention this blog prominently. So, of course, that gets embedded here.

Quantum Chemistry III

View more presentations from baoilleach.

Rotating molecules in PowerPoint

2010-05-02T09:58:00.002+02:00

Long-time MolModBasics reader/commenter NUChem posted the following question:

"One thing that I've been wondering for a while is creating movie files of optimized structures for presentations. Would it be possible for you to have a screen cast of how to take an optimized structure of cyclohexane and make a movie of the molecule rotating. Basically I have a few polycyclic molecules that I've optimized and I want to be able to have them rotating in my powerpoint presentation so my audience can see the whole molecule. Is there anyway to do this?"

I know of two ways of doing this that are relatively easy. One is to create an animated gif file using the Polyview3D web server, and inserting the file as a movie in Powerpoint. Here is a screencast of how to do this.

Another option is to set the molecule spinning in, for example, Avogadro and then simply record part of the screen using a screencasting program. Unfortunately, it looks as though you have to buy software to do this. There is a free screencasting program called Jing, but you can only save the movie file in the .swf file format, which Powerpoint can't read (and I haven't been able to find a free .swf to .mpeg converter). However, with JingPro ($15/year) you can save the file in the .mpeg format, which should work with Powerpoint, but I haven't tried it myself. I use ScreenFlow ($99) and really like it, but that only works on Macs.

UPDATE: Screencast-o-matic is free, and can make mp4 files, which can be included in Powerpoint. See this new post for more information.

Of course, with the screencast option you can add animations to Powerpoint of anything, such as vibrations, IRCs, etc. in addition to rotation. However, for more complicated animations I usually switch between Powerpoint and Jmol (more precisely a browser showing an html file with Jmol embedded) during the presentation, ever since I found out I don't have to quit the Powerpoint presentation to switch to another program.

Polarization and intermolecular interaction

2010-04-30T10:31:00.005+02:00

Figure 4.16. 0.002 au isodensity surface with superimposed molecular electrostatic potential for (a) methane using the methane density for the methane–water dimer and (b) free methane monomer. The maximum potential value, 0.01 au, is five times smaller than in previous figures to make the increase in negative charge visible. The black spheres in (a) denote the position of the three nuclei of water. The level of theory is M06/6-31+G(2d,p)//M06/6-31G(d).
Click on the picture for an interactive version.
Click here for a pop-up window
From Molecular Modeling Basics CRC Press, 2010

In a previous post I showed how the difference in the strengths of interaction between methane, methane-water, and water dimer is due to their differences in polarity, which can be visualized using molecular electrostatic potentials (MEPs). Methane interacts stronger with water than with methane because of polarization, i.e. rearrangement of the methane electrons due to the polar water molecule that, to a first approximation, induces a dipole moment in methane.

The polarization of the methane density is not visible in the MEP of the methane–water dimer (Figure 4.15b) because the polarity of the water H atom dominates the MEP in that region. But if we remove the water molecule, the net increase in negative charge in the methane molecule where it interacts with the partially positive H atom of the water is apparent (Figure 4.16).

Making the figure
This is a not a plot one makes every day, so the process is a bit involved. The main trick is to construct the methane part of the density from the corresponding localized molecular orbitals, and then tricking GAMESS into printing a file with just those LMOs that MacMolPlt can read. The procedure is described in some detail Section 5.5 of the book, so here I just post the corresponding screencast.

The book is out

2010-04-29T08:07:00.000+02:00

The book is now available for purchase, at least in the States. You can order it directly from CRC Press or, for example, from Amazon.com (though not yet Amazon.co.uk). And why just one copy? Christmas is practically around the corner.

When molecules attract

2010-04-25T09:35:00.005+02:00

Figure 4.14. M06/6-31G(d) optimized geometries of (a) methane dimer, (b) water–methane dimer where water acts as a H-bond donor, (c) water–methane dimer where methane acts as a H-bond donor, and (d) water dimer.
Click on the picture for an interactive version.
Click here for a pop-up window
From Molecular Modeling Basics CRC Press, 2010

The molecular dimers shown in Figure 4.14 have very different interaction energies: -0.5, -1.0, -0.6, and -5.1 kcal/mol, respectively; which are reasonably well reproduced at the M06/6-31+G(2d,p)// M06/6-31G(d) level of theory: -0.4, -0.5, 0.0, and -4.9 kcal/mol.

The source of this difference in intermolecular attraction can be easily visualized with electrostatic potential maps (Figure 4.15). Methane is non-polar and the main source of attraction in the methane dimer is dispersive forces (which are hard to visualize). Water is polar, and the methane–water interaction (where the water is the H-donor) is a bit stronger than the methane dimer. This is due to an electrostatic interaction - more specifically polarization, but more about this in a future post.

Figure 4.15. 0.002 au isodensity surface with superimposed molecular electrostatic potential for (a) methane dimer, (b) water–methane dimer where water acts as an H-bond donor, (c) water–methane dimer where methane acts as an H-bond donor, and (d) water dimer. The maximum potential value is 0.05 au and the level of theory is M06/6-31+G(2d,p)//M06/6-31G(d).
Click on the picture for an interactive version.
Click here for a pop-up window
From Molecular Modeling Basics CRC Press, 2010

Instructions on how to make interactive electrostatic potential maps with Jmol can be found here. Finally, I introduce a new feature (pop-up windows) to the blog because I can't figure out how to include Jmol buttons (which gives more control to the viewer) into blog posts. This feature also gives you access to the underlying GAMESS files as I have discussed here.

Double bonds - banana and otherwise

2010-04-02T13:31:00.002+02:00

Figure 4.12. (a) B3LYP/6-31G(d) optimized geometry of ethene; (b) and (c) 0.045 au isosurfaces of the two localized banana MOs corresponding to the double bond.
Click on the picture for an interactive version
From Molecular Modeling Basics CRC Press, 2010

VSEPR theory can be used to explain structures of molecules with double and triple bonds, though this is rarely done in textbooks. Figure 4.12 shows some localized MOs for ethene, where the double bond is shown to consist of two curved MOs (sometimes called banana bonds), rather than the usual sigma and pi MOs. The banana LMOs are less spread out than single-bond LMOs, leading to less repulsion and an H–C–C angle larger than 109.5°, namely 116.3° [at the B3LYP/6-31G(d) level of theory].

Figure 4.13. 0.045 au isosurfaces of a localized pi MO [(a) top view and (b) side view] and two localized sigma bond MOs in benzene. The level of theory is B3LYP/6-31G(d).

Click on the picture for an interactive version
From Molecular Modeling Basics CRC Press, 2010

In the case of benzene the LMOs actually look like sigma and pi orbitals (Figure 4.13). Figure 4.13a and b are two views of the pi-bond LMO primarily between C4 and C5. Notice, however, that there is significant delocalization onto C3 and C6. There are identical pi-bond LMOs between C3 and C2 as well as C1 and C6, and the net result is an identical C–C bond length of 1.397 Å, roughly halfway between the CC bond length in ethane (1.531 Å) and ethene (1.331 Å).

I have discussed how to make these plots in a previous post.

Canonical and localized molecular orbitals

2010-03-31T15:32:00.002+02:00

Figure 4.10. 0.045 au isosurfaces of the four valence canonical MOs of NH₃ computed using B3LYP/6-31G(d). A value of 0.045 au is chosen because it results in a 0.002 au isodensity surface when squared.
From Molecular Modeling Basics CRC Press, May 2010.
Click on the picture for an interactive version

According to Valence Shell Electron Pair Repulsion (VSEPR) theory many bond angles involving elements such as C, N, and O are close to 109.5^o because the four valence electron pairs that surround these atoms adopt a tetrahedral geometry to minimize repulsion. In the case of CH₄ the H-C-H angle is exactly 109.5^o because the repulsion between the four electron pairs in the C-H bonds are identical. NH₃ has a lone pair that is fatter than a bond near the nucleus, so the lone pair-bond repulsion is slightly larger than the bond-bond repulsion. This results in a H-N-H angle of 107^o, slightly smaller than 109.5^o.

However, this is far from obvious when looking at the four valence MOs of NH₃ (Figure 4.10) computed using B3LYP/6-31G(d). The reason is that that MOs are not unique, and that these MOs (which lead to a diagonal Fock matrix and are known as canonical MOs) are not the MOs where the electron pair repulsion is a minimum. Algorithms have been implemented that find a new set of MOs (localized MOs or LMOs) that represent a linear combination of canonical MOs for which the MO–MO repulsion is a minimum (Figure 4.11).

Figure 4.11. 0.045 au isosurfaces of the four valence localized MOs of NH₃ computed using B3LYP/6-31G(d). There are three N–H bond LMOs [(a)–(c)] and one lone pair LMO [(d)].
From Molecular Modeling Basics CRC Press, May 2010.
Click on the picture for an interactive version

LMOs for which the inter-orbital repulsion is a minimum are called energy localized orbitals or Edmiston-Ruedenberg orbitals. Other popular choices are Foster-Boys and Pipek-Mezey LMOs, which use different localization criteria.

Using MacMolPlt
The screencast below shows how to compute Ruedenberg LMOs for ammonia using GAMESS (local=ruednbrg in the $contrl group), and how to display the LMOs, as well as the canonical MOs, in MacMolPlt. I also show how to identify the HOMO and LUMO canonical MOs in MacMolPlt.

Note that I specify a geometry optimization in the GAMESS input file. GAMESS will only compute the LMOs for the optimized geometry.

Using Jmol
I use Jmol for the interactive figures and the scripts can be found here and here. Unlike the density and electrostatic potential, Jmol can generate its own grid data, so to display MO number 2 you simply use "mo 2" in the script. I use

mo 2; mo cutoff 0.045; mo fill nomesh; mo translucent 0.2

to make it a little prettier.

Jmol cannot find the MOs in a GAMESS geometry optimization file, and Jmol only stores the LMOs if present. So to display the canonical MOs with Jmol you need a single point energy calculation (output file) and to display the LMOs with Jmol you need a single point energy calculation with local=ruednberg added (output file).

Shameless book promotion

2010-03-17T18:41:00.000+01:00

Things are happening on the book front. A snappy cover (at the top of the blog) and now a promotional flyer - with table of content - which can be downloaded here. The book should be out May 10, complete with non-interactive black-and-white versions of many of the figures found in this blog.

Illustrating Entropy

2010-03-14T10:42:00.007+01:00

This screencast shows some Molecular Workbench models I made to illustrate entropy.

In the first model a monatomic gas consisting of 100 particles are initially restricted to the left container, but when the separator is removed the gas expands to fill both containers evenly, on average.

The temperature is held constant during the simulation, so the internal energy (which equal to the temperature times a constant) stays constant. The driving force behind the expansion is therefore purely entropic.

Using two simpler cases with two and three particles, I show that the probability of having N_L particles (out of N total) on the left is

where W is the weight or the number of equivalent combinations of particles in two containers. The screencast goes on to show that W is largest when the gas is evenly distributed between the two containers, so that is what happens when the separator is removed because that is the most probable thing to happen.

For N particles and two containers the number of combinations if there are N_L particles in the left container is

(as explained here and here).

As in the simpler cases, the most likely state for 100 particles is an equal distribution of particles among the two containers (N_L = 50). The main difference between 3 and 100 particles is that in the former case you'll occasionally (2/8 or 25 % of the time) see all particles in one container, whereas for 100 particles you'll never (2 x 10^-28 % of the time) see all particles in one container.

Entropy, S, is defined as

where k_B is Boltzmann's constant. This function is largest when W is largest, so one can say that the gas expands to maximize its entropy, but it's a bit like saying the gas expands because that is the most probable.

The change in entropy is

which can be approximated using the more usual equation for the entropy change upon doubling the volume of an ideal gas while keeping the temperature constant

(The latter equation becomes more accurate as the number of particles increase, because it relies on Stirling's approximation).

The model is too wide to fit in this blog but you can play around with it on this web page, or you can download the model if you have Molecular Workbench installed on your computer. Enjoy!

Related blog posts
Where does the ln come from in S = k ln(W)?
Entropy, volume, and temperature

Internal energy and molecular motion

2010-03-03T09:58:00.000+01:00

This screencast shows a Molecular Workbench model I made to illustrate the connection between internal energy and molecular motion (the model is included at the bottom of the post).

For a monatomic gas such as He or Ar the part of the internal energy (U) that depends on temperature (T) is
and has only contributions from translational motion. Kinetic theory tells us that
which, with a little algebra, gives us a direct connection between the internal energy and the average kinetic energy of the atoms in the gas:
You can see in the first container that while the internal energy reflects an average speed, the individual atoms can have very different speeds, which are constantly changing due to collisions. This can be illustrated by coloring the atoms according to the kinetic energies (red means more kinetic energy).

Increasing the temperature increases the internal energy and this is reflected in an increase in the average speed, as you can see by clicking on the "Increase T" button.

However, it is important to remember that the internal energy is a reflection of the kinetic energy, which is also a reflection of the mass. Thus, a heavier atoms will have the same internal energy and kinetic energy at the same temperature, but the atoms will move slower on average (the "Increase M" button).

Also, the internal energy can be increase by adding more particles, but this does not change the average speed of the molecules at the same temperature (the "Increase n" button).

The last ("Diatomic") button changes the atoms to diatomic molecules without changing the mass, i.e. the mass of the atoms in the second container are half that of the atoms in the first container. Because the mass is unchanged the average speed of the molecules (more precisely, of their center of mass) is the same as for the atomic gas at the same temperature, as is the translational internal energy (³/₂nRT ).

However, notice that the molecules are also rotating, meaning they have additional (internal) energy compared to the atomic gas, which amounts to nRT (³/₂nRT for non linear molecules):
The internal energy also has a vibrational contribution (here the frequency is converted to a wavenumber, i.e. units of cm^-1):
though it is hard to make out the vibrational motion in the simulation. However, for most diatomic molecules this contribution is negligible at temperatures relevant to most chemists.

You can play around with the model yourself by clicking on the picture below, or on this web page, or you can download the model if you have Molecular Workbench installed on your computer. Enjoy!

Ionic and metallic bonding

2010-02-20T21:32:00.002+01:00

Figure 4.7. 0.0005 au isodensity surface with superimposed electrostatic potential of (a) H–Li, (b) Li, (c) Li₂. The maximum potential value is 0.05 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.

Figure 4.7 highlights two other kinds of bonds. Note the decrease in the size of the Li atom (and the large positive charge) compared to the Li atom in Li–H. LiH is best thought of as Li⁺H^- (i.e., the bonding electron pair belongs to H^- rather than being shared), which is an example of ionic bonding.

Here is an interactive version of Figure 4.7a where the atomic densities of H and Li are superimposed on the density of H-Li.

Click on the picture for an interactive version

Li₂ is more polar than H₂: more positive at the ends and more negative in the middle, indicating the a larger rearrangement of electron density upon binding compared to H₂. Here is an interactive version (the corresponding Jmol script can be found here; see this post on making plots like this), where I have superimposed the electron densities of the Li atoms and Li₂:

Click on the picture for an interactive version

Looking deeper into the density (by using a larger, 0.01 au, isodensity value) reveals a density rearrangement that is quite different from H₂:

Click on the picture for an interactive version

Here you can see strong localization of electron density between the nuclei, indicating that Li₂ is best thought of as Li⁺:Li⁺, where ":" indicates an electron pair. This is an example of metallic bonding.

Both ionic and metallic bonding are distinctly different from covalent bonding, in that only covalent bonding leads to the formation of distinct molecules, while ionic and metallic bonding leads to formation of crystals.

To demonstrate this, find the lowest energy structures of H₄, Li₄, and Li₂H₂, using, for example, B3LYP/6-31G(d). Go ahead, I'll wait.

Why do I use a different isodensity value for this plot?
I usually use a 0.002 au isodensity surface, but in Figure 4.7 I use 0.0005. This is something I didn't discuss in the book but should have. When I made the plot of the 0.002 au isodensity surface with superimposed electrostatic potential of Li atom, I discovered that is was very blue, i.e. quite positive (go ahead, try it for yourself). I only got a neutral-looking atom when I went down to 0.0005 aus, which is what I have used for the figures (unless otherwise noted).

The positive electrostatic potential on the 0.002 au isodensity surface of Li indicates that a significant amount of density is found outside this surface, which is clearly not the case of the H atom, as shown in a previous post. This, in turn, indicates that the Li electron density decreases more slowly with distance compared to H, which makes sense since Li is much less electronegative than H. Thus, the 0.002 au isodensity surface is a reflection of the van der Waals surface of organic molecules (containing electronegative elements) but not alkali metals.

Covalent bonding

2010-02-18T21:14:00.001+01:00

Figure 4.5. 0.002 au isodensity surface with superimposed electrostatic potential of (a) H₂ and (b) H atom. The maximum potential value is 0.05 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.

If we look at the electron density of H₂ (Figure 4.5a), we can clearly see that at this separation the electron densities of the two H atoms have fused indicating electron sharing, a hallmark of covalent bonding. Here is an interactive version (the corresponding Jmol script can be found here; see this post on making plots like this), where I have superimposed the electron densities of the H atoms and H₂:

Click on the picture for an interactive version

What you can't see in this picture is that the electron density has rearranged significantly between the nuclei, which is the source of the bond strength. To see this we need to look at a larger isodensity value (here 0.075 au):

Click on the picture for an interactive version

One, two, three, MD

2010-02-16T08:58:00.003+01:00

Some exciting developments over at the Molecular Workbench (MW) blog run by MW author Charles Xie. I have several blog posts on MW, but it was necessary to install MW to check it out for yourself. Now it is possible to embed a MW applet in web pages (and blog posts!) like the one here (just push the play button!):

I think this is a big step forward for MW. While it is easy to download and install MW, it still removed MW a few clicks from the user and made it "appear to be yet another kind of annoying pop-up" and Charles notes.

It's very easy to do this. The screencast below shows how I made the simulation above in MW. Note that it literally takes one minute (and 5 seconds).

When you hit save you get two files: md.cml and md$0.mml. I transferred these to my web server where I had also put the MW applet (mwapplet.jar).

The html code is

<applet code="org.concord.modeler.MwApplet" archive="mwapplet.jar" height="300" width="100%"><param name="script" value="page:0:import md.cml"></applet>

To include it in a blog, where mwapplet.jar is not installed, add the server address in front of mwapplet.jar and md.cml, e.g. http://myserver.edu/md.cml.

I'm positive

2010-02-07T15:57:00.003+01:00

Figure 4.2. 0.002 au isodensity surface with superimposed electrostatic potential of (a) Li⁺, (b) Na⁺, and (c) K⁺ ion. The maximum potential value is 0.8 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.

Here is an example of how computational chemistry can be used to enhance teaching at the general chemistry level.

Why does the ionization energy decrease on from Li to Na to K? That's the same as asking why the electron affinities decrease on going from Li⁺ to Na⁺ to K⁺.

Figure 4.2 show the ions colored by how positive they are at the surface (i.e. the electrostatic potential superimposed on the 0.002 isodensity surface). The darker the
color the more positive the ion, and it is clear that the Li+ ion is “more positive” than Na⁺, which is more positive than K⁺. Thus, more energy should be released when adding an electron to Li⁺ compared to Na⁺, and hence more energy is needed to remove an electron from Li
compared to Na (and similarly for K).

The reason why Li⁺ is “more positive” than Na⁺ or, more accurately, why the potential on the 0.002 au isodensity surface is more positive for Li⁺ than for Na⁺, is that the former is a smaller ion than the latter, so the surface is closer to the +1 charge at the center of the ion.

Figure 4.4. 0.002 au isodensity surface with superimposed electrostatic potential of (a) Li⁺, (b) Ne⁺, and (c) Na⁺ ion. The maximum potential value is 0.8 au, and the level of theory is B3LYP/6-31G(d).
From Molecular Modeling Basics CRC Press, May 2010.

This rationalization is, of course, only qualitative and is is not predictive of the ionization energies among different groups. For example, based on Figure 4.4 one would expect that Ne would have roughly the same ionization potential as Na, which is not true at all.

When making these figures it is very important to get the relative sizes of the ions correct, but this can be difficult since each image can be zoomed to an arbitrary size. The screencast below shows how to control this in MacMolPlt using the Manual Windows Parameter window.

The same point can also be made with a "quick-and-dirty" electrostatic potential, as shown in this interactive figure. Here the electrostatic potential is due to a plus one charge centered at the atom (rather than the nuclear charge and the electron density) and the surfaces are the spheres defined by empirical ionic radii. The Jmol script can be found here, and the mol2 files here, here, and here.

Click on the picture for an interactive version

Tesserae, sera

2010-01-31T09:27:00.000+01:00

Figure 3.14. A solute–solvent boundary surface described by a set of interlocking atom-centered spheres. The surface is discretized by using triangular tesserae. Drawing graciously provided by Benedetta Mennucci.
From Molecular Modeling Basics CRC Press, May 2010.

Due to the computational expense of explicit solvation, implicit solvation is a popular alternative when including solvent effects in quantum chemical calculations. The polarized continuum model (PCM) and the closely related COSMO method is one of the most used implicit solvation models for QM calculations. Here the Poisson-Boltzmann equation is solved numerically by diving the solute-solvent surface into small pieces called tesserae (Figure 3.14).

Unfortunately, I have not been able to make an interactive version of this figure.

It's all about boundaries

2010-01-27T22:08:00.005+01:00

Figure 3.8. (a) Adenine (C₅H₅N₅) in a liquid drop of 246 water molecules. (b) Adenine in a periodic box of 511 water molecules.
From Molecular Modeling Basics CRC Press, May 2010.

Conceptually, the simplest way of simulating a molecule in solution is to place it in the middle of a roughly spherical ball of water molecules (Figure 3.8a) and perform an molecular dynamics simulation.

One problem is this approach is that the drop would eventually evaporate if the simulation is run long enough. Another problem with the liquid drop model is that the water molecules at the surface of the drop do not behave like water molecules in liquid water.

Therefore, most explicit solvation simulations use periodic boundary conditions (Figure 3.9).

Figure 3.9. Sketch of periodic boundary in two dimensions: (a) The position of the particles in the central box are copied and placed in neighboring boxes. Figure 3.8b shows a cube from a real simulation. (b) When a molecule tries to leave the box during an MD simulation, it reappears at the opposite end of the box, so the number of particles in the central box stays constant.
From Molecular Modeling Basics CRC Press, May 2010.

You can find interactive versions of Figures 3.8a and 3.9b here (I am grateful to Kestutis Aidas for providing the coordinates).

Click on the picture for an interactive version

And you can find an animated version of Figure 3.9 here (an example of where a movie really is worth 10,000 words).

The animation was made with Molecular Workbench (MW). You can play with the simulation here or you can download the MW file here (after you gave installed MW).

I have my moments

2010-01-23T13:28:00.003+01:00

Figure 3.7. Contour plot of the RHF/6-31G(d) electrostatic potential and 0.002 au isodensity surface of (a) CH₃COO^-, (b) HF, and (c) F₂. The maximum/minimum contour values are, respectively, 0.5/0.025; 0.1/0.005; and 0.005/0.00025 au respectively. Blue corresponds to a negative potential. In each case the outer-most contour looks like the corresponding contour in the electrostatic potential due to a charge, dipole, and quadrupole, respectively.
From Molecular Modeling Basics CRC Press, May 2010.

This figure makes 3 points:

1. It shows what the electrostatic potential of a charge, dipole, and quadrupole looks like.

2. It show the relative strengths of the electrostatic potentials due to a charge, dipole, and quadrupole.

3. It shows that the electrostatic potential of a charge, dipole, and quadrupole deviates significantly from the actual electrostatic potential near the molecular surface.

Here is a screencast showing how I made Figure 3.7a. (If I were to do it over I would have chosen red for negative and blue for positive... oh, well.)

Here is an interactive version of the figure. Click on the picture to load it. Remember it's Jmol so you can rotate it and zoom as you like (Mac users: this works best with Safari).

Figure 3.7. The RHF/6-31G(d) electrostatic potentials of acetate, HF and F₂.
Click on the picture for an interactive version

See these two posts (here and here) on how to make plots like that with Jmol. From a Jmol perspective the only new thing is that I show two isosurfaces simulateneously. This is done by loading the file twice to create two "frames" that Jmol can display simultaneously. You can find the script file with all the commands here, but the general syntax is:

load files file1.xyz file2.xyz
frame 1.1; isosurface surf1 plane {0 0 0 0} contour 40 color range -0.05 0.05 "potential.cube.gz"
frame 2.1; isosurface surf2 0.002 "density.cube.gz"

Even though you want a 2D contour plot of the potential, it is necessary to make a 3D cube file. Because I use unusually low cutoffs to show the outer contours, I had to trick MacMolPlt into making a bigger grid. I show how on the screencast below.

Common error messages in GAMESS: Failure to locate stationary point

2010-01-02T18:45:00.001+01:00

       ***** FAILURE TO LOCATE STATIONARY POINT, TOO MANY STEPS TAKEN *****
  UPDATED HESSIAN, GEOMETRY, AND VECTORS WILL BE PUNCHED FOR RESTART
**** THE GEOMETRY SEARCH IS NOT CONVERGED! ****

This error message was produced by the following input file. Can you see what's wrong?


 $contrl runtyp=optimize icharg=1 $end
 $basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

Actually, there is no problem with the input file as such. A geometry optimization is an iterative process and if the gradient it not below the convergence criteria within 20 steps, GAMESS will stop and print out the message shown above.

The solution is simply to take the last set of coordinates and run the optimization again, as I show in the screencast below.

As I've mentioned in a previous post I think the default criterion for geometry convergence (0.0001) is too strict, and the default number of steps (20) is too small. So I usually use 0.0005 and 50, respectively.

 $statpt nstep=50 opttol=0.0005 $end

Quick and dirty electrostatic potential maps

2009-12-30T21:21:00.001+01:00

In a previous post I showed how to compute an electrostatic potential map superimposed on the 0.002 isodensity surface of a molecule based on data computed using quantum chemical methods such as RHF/6-31G(d).

Previous posts (such as this one) have demonstrated that the van der Waals surface is a reasonable substitute for the 0.002 isodensity surface. Similarly, the electrostatic potential due to atomic charges can be a reasonable substitute for the electrostatic potential due to the electronic density.

The screencast above shows how to make such electrostatic potential maps using Avogadro and Jmol.

Avogadro uses an empirical method to determine the atomic charges (an integral part of the MMFF force field). It is possible to change the surface using the so-called "Iso Value" but I could not find any documentation on how that actually works. An Iso Value of 0 seems to correspond to the van der Waals sphere surface. It is currently not possible to alter the color range as far as I can see.

Jmol is not able to determine the charges, but the information can be transferred from Avogadro by saving a mol2 file. The electrostatic potential option in the Jmol menu corresponds to the following set of commands (as far as I can determine):

isosurface solvent color range -0.05 0.05 map mep
color isosurface translucent 0.5

and this set of commands can thus be used to control the color range and the nature of the surface. The default surface is a solvent accessible surface, which is slightly larger than the van der Waals surface.

A static and interactive version of the Avogadro and Jmol electrostatic potential map, respectively, can be found here.

Electrostatic potential map made with Avogadro
Click on the picture for an interactive version made with Jmol

Steric strain vs electrostatic attraction

2009-12-29T15:37:00.001+01:00

Figure 3.5. (a) RHF/6-31G(d) 0.002 au isodensity surface with superimposed electrostatic potential for (a) cis-HO(H)C=C(H)OH and (b) cis-CH3(H)C=C(H)CH3 and. In both cases, the maximum potential value is 0.05 au. (click on the picture for a bigger version).
From Molecular Modeling Basics CRC Press, May 2010.

This figure shows how the electrostatic potential superimposed on the 0.002 au isodensity surface can be used to rationalize why cis-HO(H)C=C(H)OH is more stable than trans-HO(H)C=C(H)OH, while the opposite is true for CH3(H)C=C(H)CH3.

The figure clearly shows the difference in polarity between hydroxyl and methyl groups. For the OH substitutent case the positive (O)H atom is close to the negative O(H) atom in the cis isomer. The CH3 group is non-polar and larger and the overlapping density indicates steric strain.

The figure is made with MacMolPlt as described in a previous post. Below is an interactive version made with Jmol (as described in a previous post). You'll notice that the color scheme is somewhat different, because Jmol maps the electrostatic potential value to a color in a different way than MacMolPlt. However, the general conclusion drawn from both programs is clearly the same.

Figure 3.5. (a) RHF/6-31G(d) 0.002 au isodensity surface with superimposed electrostatic potential for (a) cis-HO(H)C=C(H)OH and (b) cis-CH3(H)C=C(H)CH3 and.
Click on the picture for an interactive version

Electrostatic potential maps reloaded

2009-12-27T09:58:00.003+01:00

Here is an interactive version of a figure I described in a previous post. Click on the picture to load it. Remember it's Jmol so you can rotate it and zoom as you like (Mac users: this works best with Safari).

Figure 3.4. The RHF/6-31G(d) electrostatic potential of water.
Click on the picture for an interactive version

The Jmol animation loads a file (fig3-2-2.xyz), which you can access by right-clicking on the Jmol animation once you have loaded it (as described in a previous post). This file also contains all the necessary commands.

Figure 3.4.c is the most common depiction of electrostatic potential maps and the Jmol general syntax for the command for this is

isosurface 0.002 "density.cube.gz" color range -0.05 0.05 "potential.cube.gz"

The screencast below shows how I made the cube files that contain the electron density and electrostatic potential information using MacMolPlt. The first part is identical to the screencast in a previous post on obtaining the electron density cube file.

The program I used to convert the MacMolPlt file to a cube was written by Jonathan Gutow and can be download here.

I start by loading the RHF/6-31G(d) optimized geometry I computed for a previous post, and reorienting. The orientation makes it easier to define the plotting plane for the contour plots.

Wishlist

2009-12-26T10:25:00.003+01:00

Christmas is only 363 days away. It's never to early to start working on a new wish list.

[source]

Electrostatic potential maps

2009-12-21T09:35:00.003+01:00

Figure 3.4. (a) Contour plot of the electrostatic potential of H₂O. The maximum and minimum contour value is 0.5 Hartrees/electronic charge (au). (b) The corresponding 0.05 au isopotential surface. (c) The electrostatic potential displayed on the 0.002 au isodensity surface of water. The maximum (darkest blue) value corresponds to 0.05 au.
From Molecular Modeling Basics CRC Press, May 2010.

Here is a screencast of how I made the figure:

The files I used were created in a previous post. The various values specified above are determined using trial and error, i.e. I kept fiddling with it until it looked good to me. When using plots like this it is important to specify these values, because using different values can lead to very different looking plots.

Also, different programs use different color scales and color intensity to indicate positive and negative charge, so it is rarely possible to directly compare electrostatic potential maps from different programs.

See this post about using screen capture to get a file with the graphic. This was done for each plot, and the combined in a word processor.

The interactive version of this figure is the subject of a future post.

Electron density doesn't always tell the whole story

2009-12-19T11:00:00.002+01:00

Figure 3.3. (a) RHF/6-31G(d) 0.002 au isodensity surface and (b) van der Waals surface of cis-HO(H)C=C(H)OH (click on the picture for a bigger version).
From Molecular Modeling Basics CRC Press, May 2010.

Looking at the electron density or the van der Waals surfaces one would expect steric strain in cis-HO(H)C=C(H)OH. (See this post on how to make such pictures and this post on how to make an interactive version such as the one shown below.)

However, contrary to CH₃(H)C=C(H)CH₃, the cis isomer is more stable than the trans isomer. This is of course due to the electrostatic attraction of the O and H atom in the hydrogen bond. Here the electron density tells only part of the story because this attraction involves the nuclei as well. Here an electrostatic potential map is more useful. How to make such plots will be the subject of future posts.

Figure 3.3. The RHF/6-31G(d) electron density of cis-HO(H)C=C(H)OH.
Click on the picture for an interactive version

Common error messages in GAMESS: charge and multiplicity

2009-12-11T14:05:00.006+01:00

  *** CHECK YOUR INPUT CHARGE AND MULTIPLICITY ***
THERE ARE    33 ELECTRONS, WITH CHARGE ICHARG=  0
BUT YOU SELECTED MULTIPLICITY MULT=  1

This error message was produced by the following input file. Can you see what's wrong?


 $contrl runtyp=optimize $end
 $basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

As GAMESS tells you, it is not possible to have singlet state (i.e. a multiplicity of 1) if you have an odd number of electrons. You need to tell GAMESS whether you want to change the multiplicity (for example, $contrl mult=2 scftyp=rohf $end) or add ($contrl icharg=-1 $end) or remove ($contrl icharg=1 $end).

When looking at the molecule it is pretty clear that you want to remove an electron, since ammonium has a positive charge. So the first line in the input should be changed to

 $contrl runtyp=optimize icharg=1 $end

Click on the picture for an interactive version

Notice that you don't have to specify that the positive charge is on ammonium. The quantum mechanics will take care of that.

Electron density and steric strain

2009-12-09T13:58:00.004+01:00

Figure 3.2. (a) RHF/6-31G(d) 0.002 au isodensity surface and (b) van der Waals surface of cis-CH₃(H)C=C(H)CH₃ (click on the picture for a bigger version).
From Molecular Modeling Basics CRC Press, May 2010.

This figure shows how the electron density and van der Waals surfaces can be used to visualize steric strain. See this post on how to make such pictures and this post on how to make an interactive version such as the one shown below.

Figure 3.2. The RHF/6-31G(d) electron density of cis-CH₃(H)C=C(H)CH₃.
Click on the picture for an interactive version

Common error messages in GAMESS: error reading IDUM

2009-12-08T21:13:00.001+01:00

 **** ERROR READING VARIABLE IDUM     CHECK COLUMN  1
H     1.0    -0.08386     2.76975     0.70945                                
....V....1....V....2....V....3....V....4....V....5....V....6....V....7....V....8

This error message was produced by the following input file. Can you see what's wrong?


 $contrl runtyp=optimize $end
$basis gbasis=pm3 $end
 $data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

The answer is the missing space in front of the $basis, and the solution is to add a space in front of it, just like for $contrl and $basis.

Why does GAMESS complain about a line in $data when the problem is with $basis? The reason is that it is possible to specify a basis set for each atom in $data, so GAMESS tries to read this information if no $basis group is found.

The error message would clearly benefit from the line: "Did you mean to specify a $basis group?"

The future of scientific publishing is here today

2009-12-06T11:00:00.000+01:00

I recently came across this paper on embedding interactive 3D graphics in pdf files. Unfortunately, I couldn't animate the 3D figures, and I wrote the author (Vlad Vasilyev) who kindly wrote back that:

(1) the publisher (Springer) is a bit behind the times so the interactive figures are found in the supplemental material, but

(2) the ACS is more on the ball, as can be seen in this paper (this is the one I use in the screencast), and

(3) he has made a bunch of tutorials on how to create such diagrams on this page.

This page also contains some pdf files with 3D graphics you can play around with if you don't have access to the papers. Be sure to use Adobe Acrobat Reader 9.

Note that it is also possible to include animations such as vibrational motion (NB: 16 MB file)

Very cool! I haven't tried this myself, but you can be sure I will. Unfortunately, it seems that you need a piece of commercial software (Adobe Acrobat 9 Pro Extended) to do this.

On a related note, check you the interactive figure of interactive figures created by fellow blogger Henry Rzepa (please give the page some time to load).

Electron density reloaded

2009-11-29T15:26:00.004+01:00

Figure 3.1. The RHF/6-31G(d) electron density of water.
Click on the picture for an interactive version

The Jmol animation loads this file (h2oprinc.xyz), which looks like this

3
jmolscript: script "http://propka.ki.ku.dk/~jhjensen/h2odensity.spt"
O -0.0000 0.0643 0.0000
H -0.7541 -0.5091 0.0000
H 0.7541 -0.5091 0.0000

and which, in turn, loads a script (h2odensity.spt), which looks like this

isosurface planex plane {0 0 0 0} contour 20 color absolute 0.002 0.05 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
delay 3
spin y 20
delay 10
spin off
isosurface planey plane {1 0 0 0} contour 20 color absolute 0.002 0.05 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
spin y 20
delay 10
spin off
isosurface planez plane {0 1 0 0} contour 20 color absolute 0.002 0.05 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
spin x 40
delay 10
spin off
isosurface threed 0.002 "http://propka.ki.ku.dk/~jhjensen/h2oprinc.cube.gz"
color isosurface red ; color isosurface translucent 0.15
spin y 20
delay 10
spin off
color isosurface red ; color isosurface translucent 0.5
select all; spacefill 100 %babel
spin y 20
delay 10
spin off

The screencast below shows how I made the cube file that contains the electron density information using MacMolPlt. The program I used to convert the MacMolPlt file to a cube was written by Jonathan Gutow and can be download here. In the screencast I mistakenly named the cube file h2odensity.cube.gz, but that's easy to change.

I start by loading the RHF/6-31G(d) optimized geometry I computed for a previous post, and reorienting. The orientation makes it easier to define the plotting plane for the contour plots.

Common error messages in GAMESS: No $data

2009-11-25T08:36:00.000+01:00

**** ERROR, NO  $DATA   GROUP WAS FOUND
EXECUTION OF GAMESS TERMINATED -ABNORMALLY-

This error message was produced by the following input file. Can you see what's wrong?


 $contrl runtyp=optimize $end
 $basis gbasis=pm3 $end
$data
Title
C1
N     7.0    -0.39094     1.95659     0.14008
H     1.0     0.38874     1.60529    -0.40413
H     1.0    -0.08386     2.76975     0.70945
H     1.0    -0.72485     1.22934     0.80007
H     1.0    -1.14035     2.30329    -0.48754
O     8.0    -0.64579     0.16732     2.03360
H     1.0    -0.26212    -0.73042     2.10569
H     1.0    -1.00756     0.26750     2.93979
O     8.0    -1.80535     3.31298    -1.59619
H     1.0    -1.39440     3.81065    -2.33214
H     1.0    -2.74148     3.57968    -1.71559
O     8.0     0.26578     4.05264     1.54485
H     1.0     1.03270     4.27032     2.11226
H     1.0    -0.26135     4.87344     1.64760
 $end

The answer is the missing space in front of the $data, and the solution is to add a space in front of it, just like for $contrl and $basis.

Electron density

2009-11-21T13:00:00.005+01:00

Figure 3.1. (a) A contour plot of the density of H₂O computed using RHF/6-31G(d). The maximum and minimum contour values are 0.05 and 0.002 aus. (b) The corresponding 0.002 au isodensity surface. (c) The surface corresponding defined by atomic spheres with van der Waals radii.
From Molecular Modeling Basics CRC Press, May 2010.

Here's a screencast about how I made the figure:

When making the contour plot (Figure 3.1.a) I pick 0.05 au as the maximum value (this will be the contour line closest to the nuclei) and 25 contour lines. This means that the spacing between the contour lines and thus the outer contour line will be 0.05/25 = 0.002 au.

Another, more common, representation of the density is a 3D version of one of the contour values: the isodensity surface (Figure 3.1.b). A common choice is 0.002 au, since that corresponds roughly to experimental estimates of molecular size, such as the van der Waals surface (Figure 3.1.c).

Unfortunately, MacMolPlt doesn't have a van der Waals display style, so I have to use Avogadro. This means I have to re-size (by eyeball) this part of the figure to make it the same size as the density plots. See this post about using screen capture to get a file with the graphic.

The interactive version of this figure is the subject of a future post.

The book and color figures

2009-11-21T12:18:00.000+01:00

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

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

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

A sense of scale

2009-11-11T21:00:00.000+01:00

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

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

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

The trouble with transition metals II: SCF convergence

2009-11-09T22:03:00.003+01:00

In a previous post I showed how to build two structural models containing transition metals: a small model of a zinc finger and a ferrocene. This post is about computing the energy using GAMESS.

The first challenge in computing the energy is choosing the correct charge and multiplicity. Most organic molecules are singlets (i.e. all orbitals are doubly occupied) and have an easily identifiable charge. This is often not the case with molecules containing transition metals. Many transition metal atoms have unpaired d electrons and can have several different charges (oxidation states). Also, the charge of the ligands is not always obvious, though they almost always are singlets.

The zinc finger model is relatively easy: zinc is almost always Zn²⁺ and has 5 doubly occupied d orbitals, so that's a pretty safe bet. The zinc finger model contains two neutral groups (the imidazole rings) and two negatively charged groups (CH3-S^-). So the charge of the entire system is zero and the multiplicity is 1.

The ferrocene model is a bit more complicated: iron is commonly found in both the +2 and +3 oxidation state. Before one starts the calculation it is important to find out which one is appropriate for ferrocene. Google says +2 (meaning 6 d electrons) and all doubly occupied orbitals (see the orbitals at the bottom of the page), so a singlet. Furthermore, it is also important to know that each ring has a -1 charge, so the total charge of the complex is 0 and the mutiplicity is 1.

Now, I was all set to talk about SCF convergence problems, i.e. that the zinc finger would converge fine but the ferrocene would give problems. But when I ran a RHF/6-31G(d) energy calculation

$contrl runtyp=energy $end
$basis gbasis=n31 ngauss=6 ndfunc=1 $end
$scf dirscf=.t. $end

both runs converged fine!

"Luckily" when using a smaller basis set (RHF/STO-3G) leads to convergence problems - for both. Here is the relevant part of the SCF output for zinc finger.

                                                                                                                                                                     NONZERO     BLOCKS
ITER EX DEM     TOTAL ENERGY        E CHANGE  DENSITY CHANGE     ORB. GRAD      INTEGRALS    SKIPPED
 1  0  0    -3064.6703808712 -3064.6703808712   1.733010527   0.000000000        9751402     457807
 2  1  0    -3059.0602240553     5.6101568159  16.724847696   0.794728909        9737106     471195
 3  2  0    -3053.0258026844     6.0344213709  16.846293078   0.656307634        9778274     459644
 4  3  0    -2935.9005712044   117.1252314800 115.104028436   1.822663152        9792060     446482
 5  4  0    -2923.8944526114    12.0061185930 115.740577390   0.620923864        9792989     441267
 6  5  0    -2768.2361993103   155.6582533011 115.747362137   5.442595435        9783876     443565
 7  6  0    -2921.1390148151  -152.9028155048 115.772970667   0.810505145        9782321     444765
 8  7  0    -2775.0326048138   146.1064100013 115.770983819   4.741230605        9783214     444582
 9  0  0    -2920.9171138162  -145.8845090024  60.570511833   0.796143103        9783475     444487
10  1  0    -3003.9501726068   -83.0330587906  54.636432025   0.749401761        9790866     444493
11  2  0    -2862.1853212282   141.7648513786 115.095938299   1.824897836        9796139     439522
12  3  0    -2921.4368542174   -59.2515329892 115.495602218   0.722894426        9783664     443903
13  4  0    -2848.2964423996    73.1404118178 115.538938253   2.394847457        9782038     444795
14  5  0    -2905.4155965583   -57.1191541586 115.529755578   0.678236838        9781255     445351
15  0  0    -2859.1068450197    46.3087515386  14.660262004   1.964358011        9783437     444643
16  1  0    -3019.7713654862  -160.6645204665  14.032213984   0.801579670        9794434     443563
17  2  0    -2934.9445251088    84.8268403773  99.891887745   1.748565491        9794346     445610
18  3  0    -2905.6599083961    29.2846167127 100.392269592   1.167341691        9789953     445157
19  4  0    -2864.9366576173    40.7232507788 104.001625891   1.797814244        9781805     446591
20  5  0    -2855.1081148594     9.8285427580 103.862237046   1.447532031        9781260     446026
21  6  0    -2865.3401581051   -10.2320432457  82.200569407   1.143844678        9780463     446052
22  0  0    -2857.7261418857     7.6140162194  78.128990865   1.766124007        9781544     445447
23  1  0    -3036.7926259046  -179.0664840190   5.185277571   0.905312779        9796083     442344
24  2  0    -2977.1610659805    59.6315599241   5.419478057   1.193389134        9796601     444172
25  3  0    -2859.5939849458   117.5670810347 114.799971576   1.724872824        9789000     444821
26  4  0    -2913.7281172708   -54.1341323250 115.553208606   1.094842493        9785096     443663
27  0  0    -2845.9860251842    67.7420920865  16.890315601   1.995942532        9783715     444171
28  1  0    -2989.2091196251  -143.2230944409  16.260193167   1.244487455        9796844     439846
29  2  0    -2911.6240721064    77.5850475187 114.870828357   1.711669748        9798193     439590
30  3  0    -2897.8498052725    13.7742668339 115.443658127   0.704588469        9789911     442464

SCF IS UNCONVERGED, TOO MANY ITERATIONS
  TIME TO FORM FOCK OPERATORS=      84.7 SECONDS (       2.8 SEC/ITER)
  FOCK TIME ON FIRST ITERATION=       2.6, LAST ITERATION=       3.0
  TIME TO SOLVE SCF EQUATIONS=       0.4 SECONDS (       0.0 SEC/ITER)

FINAL RHF ENERGY IS        0.0000000000 AFTER  30 ITERATIONS

Clearly, the this SCF is not converging and giving it more iterations is not going to make the problem go away. Instead, different algorithms for SCF convergence are needed and

$scf dirscf=.t. diis=.t. damp=.t. $end

leads to convergence in for both ferrocene and zinc finger.

DIIS stands to Direct Inversion of the Iterative Subspace, which is an extrapolation procedure, while DAMP "damps" any oscillations between in the SCF.

There are a few other tricks to SCF convergence but they are clearly not needed here. If DIIS and DAMP doesn't solve the problem for you, leave a message describing your molecule and I'll see what I can do.

The trouble with transition metals I: molecule building

2009-11-01T11:58:00.005+01:00

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

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

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

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

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

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

Get a reaction: Intrinsic Reaction Coordinate

2009-10-25T16:54:00.000+01:00

Lucas requested a post on computing an intrinsic reaction coordinate (IRC). An IRC is a special case of the minimum energy path (MEP) connecting reactant, transition state, and product (more about this at the end of the post).

There are many uses for the IRC, the most basic of which is to verify that the transition state you found actually connects the reactants and products you think it does (this is not always obvious from the normal mode associated with the imaginary frequency). It is also a good way to identify any reaction intermediates you haven't thought about. The IRC is also a way generate an animation of the reaction (see for example ChemTube3D).

An IRC is generated by following the gradient down-hill from the transition state to products and reactants. (It's a bit like a geometry optimization initiated at the transition state, but the algorithm used to follow the gradient is different). Computing a complete IRC thus requires a minimum of two separate runs, and you must find the TS before you can compute an IRC.

The screencast shows how to generate the IRC associated with the 3TSa transition state I found in a previous post. The first input file I generated contain the 3TSa coordinates and the following keywords:

$contrl runtyp=irc $end
$basis gbasis=pm3 $end
$irc saddle=.t. tsengy=.t. forwrd=.t. npoint=100 stride=0.2 opttol=0.0005 $end

saddle=.t. tells GAMESS that the coordinates in $data corresponds to a TS. Because the TS is a stationary point the gradient is zero (or very small), so there is no gradient to follow down-hill. Therefore, the IRC is initiated by displacing the geometry along the normal mode corresponding to the imaginary frequency. Therefore, saddle=.t. GAMESS expects a $HESS group, which contains the frequency information.

Whether the normal mode points towards the reactants or products is completely arbitrary, and the forwrd=.f. keyword is used to change the direction in which the normal mode points. forwrd=.t. is the default and I only include it in the first run for completeness.

tsengy=.t. tells GAMESS to re-compute the energy of the TS (re-compute because you already computed it when you computed the frequencies). If you want to make an energy plot in MacMolPlt as I show in the screen cast, you must set tsengy=.t.

The displaced geometry will have a non-zero gradient, the negative of which points downhill. The geometry is displaced in this direction, the gradient is recomputed, and the procedure is repeated npoint=100 times or until the (root-mean-square) of the gradient is below opttol=0.0005 atomic units, at which point you are at a minimum. (Notice that the default for npoint is 1, and must always be changed).

stride=0.2 tells GAMESS how far along the gradient it should displace at each step (in this case 0.2 atomic units). The default is 0.3, but when I tried this, the IRC calculation in the forward direction failed after 7 steps, because the coordinate deviated too much from the minimum energy path, so I had to decrease it. This means I will have to take more steps to reach the minimum.

As you saw in the screencast, GAMESS takes 100 steps in each direction, which means the IRC does not go all the way to the minimum on each side of the TS. GAMESS does print out restart information to continue the IRC, which I didn't do in this case.

Finally a note on MacMolPlt. As I show in the screencast, MacMolPlt can be used to combine the two IRC runs (forwrd=.t. and forwrd=.f.) by using Files > Add Frames from File... You should first load the file that goes to product and then add the file from the IRC that goes to reactants. Here you should tell MacMolPlt to reverse the order of the structures in the file, so the animation goes from reactions to transition state to product.

This is done by selecting "make these points negative". "Negative" refers to the x-coordinate in the energy plot. The x-coordinate is what separates an IRC from a minimum energy path. An IRC is a plot of the energy of each structure along the minimum energy path vs the root-mean-square change in the mass weighted Cartesian coordinates relative to the TS structure. The TS thus has a x-coordinate of zero, while structures leading to the reactants are defined by negative RMSD values.

If the energy is plotted against something else, such as a bond-length, then it is simply a minimum energy path. Because an IRC uses mass-weighted coordinates, using different isotopes will lead to different IRCs, and one use of IRC is to determine tunneling corrections to the reaction rate.

Avogadro 1.0 Released

2009-10-24T17:52:00.004+02:00

Avogadro - Code Swarm from Marcus Hanwell on Vimeo.

Avogadro 1.0, i.e. the first non-beta version, was released yesterday. The best place to read about this is on Marcus Hanwell's blog (where the movie comes from) and Geoff Hutchison's interview on SourceForge.

The movie (which is best viewed on the vimeo site using full screen) very elegantly shows the evolution of the program (notice the time on the lower right hand side of the screen), from the very beginnings in 2006 when Geoff and Donald Curtis started the project.

For the historians in the audience, Donald has previously worked with me on Ghemical-gms. Ghemical is an open-source model builder based on Babel written by Tommi Hassinen, which Donald made installable on Windows and Mac's to which he added an interface to GAMESS. During the process Donald and Geoff hooked up and the rest, as they say, is history.

Proteopedia

2009-10-18T11:39:00.000+02:00

Proteopedia is just what the name suggests: a wikipedia for proteins. What makes it specific to proteins is the integration of Jmol for animated and interactive figures. It has a long list of pre-made entries ready for use in a lecture and, like a real wiki, you can add content as well.

A lot of the features are explained using a set of very nice screencast. The screencast above is simply a recording of the first 2 minutes of the first video (protopedia does not offer a way to embed the videos) narrated by Eran Hodis.

Ribosome Animation

2009-10-14T19:14:00.001+02:00

For a while I thought I would be only science blog that would make it through October without mentioning the Nobel Prize. But via Chemical Blogspace I came across this video interview of the 2009 Nobel Prize co-winner in Chemistry Venki Ramakrishnan posted on The Sceptical Chymist (note: spell-check when you name your blog!):

During the interview they show some stunning molecular animations of the ribosome in action, which are credited to Said Sannuga and Ramakrishnan. I wonder ... Sure enough there are several very interesting animations on Ramakrishnan's web site. The screencast above shows a few seconds of one of them. For the rest, go to the site and have a look.

Building a complicated molecule: meet Marvin

2009-10-04T12:52:00.001+02:00

In a previous post I showed how to use PubChem Editor to draw a complicated molecule and import it into Avogadro using the insert SMILES option. Unfortunately, Avogadro does not interpret the chirality information correctly, and chiral structures often have to be fixed.

In the screencast above I show another program, Marvin Sketch, that allows you to draw a structure in 2D and convert it (correctly) to 3D. The coordinates can then be transferred to Avogadro by saving a .mol file.

I also use Marvin to draw 2D structures for the blog. For example

However, Marvin Sketch is much more than a drawing program. The tools menu contains many powerful cheminformatics tools such as pKa prediction, isomer search, etc.

Marvin Sketch can also be download as a stand-alone program free of charge for academics.

Look who's talking...the ACS

2009-10-02T08:10:00.004+02:00

The ACS has released audio recordings synchronized with powerpoint slides of some of the presentations given at the August ACS meeting. A list of the presentations can be found here.

So far, it appears to be free of charge. Let's hope they keep it that way.

Two sets of talks are related to this blog: molecular visualization and animation in chemical education, and I list them here.

I was particularly interested to learn that a GAMESS/VMD interface is well on its way (talk 130-3).

130. Molecular Visualization
130-1 Dr. Michelle M. Kuttel: Visualization of cyclic and multibranched molecules with VMD

130-2 Dr. Katrin Stierand: PoseView: 2D Visualization of protein-ligand complexes

130-3 Dr. Jan Saam: A general interface to quantum chemistry simulations in VMD

130-4 Dr. Randall B. Shirts: Boltzmann 3D simulations for visualizing molecular motion in the classroom and laboratory

130-5 Dr. Ugo Varetto: Molekel: A program for the visualization of quantum chemistry data

130-6 Dr. William F. Polik: WebMO: Web-based, state-of-the-art, and cost effective computational chemistry

130-7 Dr. Maciej Haranczyk : Visualization of Molecular Orbitals and the Related Electron Densities

340-2 Dr. Michel F. Sanner: Molecular visualization and animations using PMV

340-3 Dr. Thomas Ferrin: UCSF Chimera

180 Using Technology to Enhance Learning in Organic Chemistry

180-1 Dr. Matthew R. Radcliff: Using animation to improve visualization and understanding of reaction mechanisms

180-2 Dr. Ghislain Deslongchamps: Organic Chemistry Flashware: Animating mechanisms and orbital interactions

180-3 Dr. Philip A Janowicz: Creation and implementation of a fully online organic chemistry course

180-4 Dr. Suzanne M. Ruder: Using technology to engage the organic chemistry student