Melting: a simple model

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

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

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

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

Many-body effects in the water trimer

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

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

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.