Amide hydrolysis, revisited 2

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

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

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.

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.

It takes a village to solve a Jmol puzzle

(Click-drag on the image to rotate the molecule!)

The past few days I have been trying to put Jmol animation into the blog. I knew this was possible because of the pioneering efforts of Felix over at Chemical Quantum Images.

His post on vibrations in p-cresole has been a long-time inspiration. I found the entry through google because I was interested in animating vibrations with Jmol, and it taught me a lot about that.

But more importantly, it was the first blog I came across where someone talked about the details of their computations and showed images resulting from the computations simply because they were beautiful. So Chemical Quantum Images was a major inspiration for this blog. Thanks, Felix!

And now, through the comments to that post, and some very helpful suggestions by Noel over at Noel O'Blog (thanks, Noel!), I now also know how to embed Jmol within this blog. Noel pointed me to another way of doing it, via a Microjmol Widget application written by Jeffrey Moore and that provided me with the penultimate piece of the puzzle.

Anyway, just in case anyone else is interested in doing this: the text I put into this blog window to generate the Jmol window can be found below. Note that this will not show up in Preview, only when you post (and that was the final piece of he puzzle).

<script src="http://propka.ki.ku.dk/~jhjensen/Jmol.js" type="text/javascript"></script> <script type="text/javascript"> jmolInitialize("http://propka.ki.ku.dk/~jhjensen"); jmolApplet(300, "load http://propka.ki.ku.dk/~jhjensen/cyclohexane_movie.xyz;");jmolBr();jmolButton("spin y" ,"spin y");jmolButton("spin off" ,"spin off");<br /></script>

Some Jmol basics

Here is a Jmol application I wrote to illustrate rotational and vibrational energy states (in HCl) during a p-chem level lecture on energy states and statistical mechanics. You can find it here.

There are 3 main points:

1. It's an example of how you can use Jmol to visualize molecular motion. Maybe you like it and want to use it, or get inspired to do something much cooler.

2. The vibrational mode comes from a GAMESS calculation, and I show how you can access the GAMESS output file directly from Jmol. (Note you can get the coordinates of any molecule displayed with Jmol that you find on the web that way.)
November, 2010: In newer versions of Jmol the menu has changed a bit.  You now access the coordinates from "About" at the very bottom of the menu.

3. As with many web pages, the underlying html can be accessed from the browser.

So you have everything you need to reconstruct this example if you want. Of course you need to install Jmol on the computer that hosts the web site, and the GAMESS file needs be in the same directory as Jmol (at least on my web server).

The increases in rotational speed (100 and 1000) and vibrational displacements (0.1, 0.5, and 1.0) that I chose in the html code are completely arbitrary. I thought this got the point across. Click here for a complete list of commands.

A more complicated example (H2O) can be found here.