Overview

The binding energy (E_bind) of two interacting molecules can be determined using the following equation:

\[ E_{\mathrm{bind}} = E_{\mathrm{dimer}} - \left ( 2 \times E_{\mathrm{monomer}} \right ) \]

where E_dimer is the energy of the dimer and E_monomer is the energy of the monomer. Because two water molecules favorably interact, there is a lowering of the potential energy. By subtracting the energy of two, non-interacting optimized water molecules from the energy of optimized dimer system, this gives the energy of the interaction (i.e. the binding energy).

Therefore, we must optmize a single, isolated water molecule as well as the dimer complex. We will examine the low-energy non-planar open C_s conformer of the water dimer (pictured below).

Setup and Run

Create the following directories.

• ~/chem/water/monomer/b3lyp/dz/
• ~/chem/water/dimer/b3lyp/dz

Within each directory, create input files, one for the water monomer and one for the water dimer. Be sure to include a copy of the PBS file in each directory.

• Water Monomer Input File
• Water Dimer Input File
• PBS File

You will notice that the molecule specification in the input file is given as a Z-matrix and not as Cartesian coordinates. The Z-matrix is carefully crafted to ensure the C_2v symmetry of the monomer and the C_s symmetry of the dimer.

Submit the job to the queue.

% qsub g09.pbs

Your job should now be in the queue. To check if it has been queued up, type

% qstat | grep $USER

This will only show your jobs that are in the queue. Once the job completes, it will no longer be in the queue. A 1.out file should be in your directory if the job ran.

Analyze the Output File

Water Monomer

The output file contains a plethora of information. For now, we want to pull the final electronic energy from the output file.

% grep "SCF Done:" 1.out

The last energy is what we want. Another command to type to simply get the last energy is

% grep "SCF Done:" 1.out | tail -1
 SCF Done:  E(RB3LYP) =  -76.4206270785     A.U. after    1 cycles

You must also check the nature of the stationary point from the frequency computation. The optimized geometry should be a minimum on the potential energy surface. Verify this by ensuring all the vibrational modes are real (not printed with a negative sign).

% grep "Frequencies" 1.out
 Frequencies --   1658.8679              3751.3292              3852.7337

Here, all the frequencies are real (reported as positive numbers). Therefore, the geometry is indeed a minimum (on the B3LYP/cc-pVDZ potential energy surface).

Water Dimer

Do the same with your water dimer output file.

% grep "SCF Done:" 1.out | tail -1
  SCF Done:  E(RB3LYP) =  -152.854436576     A.U. after    1 cycles

We now have our energies to calculate the binding energy.

Binding Energy

Let us convert our absolute electronic energies into binding energies according to the equation at the beginning of this tutorial. I will then convert this energy to kcal mol^–1 by multiplying by 627.51. This can be done in an Excel spreadsheet for convenience.

Here we see that the binding energy for our water dimer is –8.27 kcal mol^–1 on the B3LYP/cc-pVDZ potential energy surface. The binding energy is negative indicating an attractive interaction.

Other Binding Energies

The experimental value for the binding energy of the water dimer is –5.4±0.7 kcal mol^–1.(Curtiss, Frurip, and Blander 1979) Below are binding energies of the water dimer as fully characterized at the listed methods and basis sets.

^aNumber of basis functions for the water dimer (divide by 2 for water monomer)

Below is a plot of the data in the table above. The solid black line represents the experimental value and the dashed black lines represent the experimental error.

Questions

• Which levels of theory agree well with experiment?
  
• What levels of theory would you not recommending using?
• What other factors should you consider when comparing these results to experiment?
  
• Would you expect these computations to agree well with experiment? How can one make better predictions?