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.

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.

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.

A typical set of GAMESS calculations

In this post I show the sequence of GAMESS calculations necessary to compute the free energy of the water dimer molecule at the B3LYP/6-31G(d)//PM3 level of theory. This notation means that the geometry and the free energy contribution is obtained with PM3, while the electronic energy is computed using B3LYP/6-31G(d) and the PM3 geometry. (By default the free energy is computed at 298.15 K and 1 atmosphere pressure.)

The steps are as follows:

1. Build the molecule

2. Optimize the geometry using PM3

3. Compute the frequencies for the optimized geometry to
a. verify that you found a minimum
b. compute the sum of the translational, rotational, and vibrational free energies at the PM3 level.

4. Compute the B3LYP/6-31G(d) energy for the PM3 optimized geometry

Some of the keywords I used have been described in a previous post, and I actually copy the keywords from this post to save typing.

The electronic energy for the water dimer is -152.7533 atomic units (au) while the sum of the translational, rotational, and vibrational free energies is 13.820 kcal/mol. If you repeat these calculations for the water molecule you get -76.3715 au and 2.271 kcal/mol (1 au = 627.51 kcal/mol).

So, the change in electronic energy on going from two water molecules to the water dimer is

ΔE = (-152.7533 - 2(-76.3715))*627.51 = -6.4 kcal/mol

The corresponding free energy change is

ΔG(298K) = -6.4 + 13.820 - 2(2.271) = -6.4 + 9.3 = 2.9 kcal/mol

Here I make the usual assumption that the electronic free energy is the electronic ground state energy.

Notice that the free energy change is positive meaning the water dimer is not predicted to form at this temperature and pressure (1 bar) under equilibrium conditions (the relative probability of observing the water dimer can be computed with this equation). This is because of the entropy of loss on going from 2 particles to 1.

The AutoOpt tool in Avogadro

I am have been meaning to do a screencast of the Auto Optimization feature in Avogadro for a while. Fortunately Geoff Hutchison beat me to it, as you can see above, and did a much better job than I could have. To prove that point I made the screencast below. It shows how the autoopt tool can be used to illustrate intermolecular interactions - in this case hydrogen bonding between water molecules. Note how the H atoms on one water follow the O atom on the other when you move it (electrostatic attraction), how it's impossible to get the O atoms next to each other (electrostatic repulsion), and how one molecule gets out of the way when you try to push the other too close (steric repulsion).