Polarization and intermolecular interaction

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

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

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

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

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

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

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