# Water Dimer

Characterize the binding energy of the water dimer

This exercise will show you how to determine the binding energy for the water dimer in the non-planar open Cs conformation .

## Overview

The binding energy (Ebind) 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 Edimer is the energy of the dimer and Emonomer 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 Cs 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.

WATCH: Setup the directory and input file for the water monomer.

Commands:

• mkdir -p ~/chem/water/monomer/b3lyp/dz to make the directory from your HOME directory
• cd ~/chem/water/monomer/b3lyp/dz to change into the directory
• vi 1.inp to create the input file
• i to enter Insert Mode
• right click (or middle click) to paste
• Enter to insert a blank line at the end of the file
• ESC to enter Command Mode
• :wq to save and close the 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 C2v symmetry of the monomer and the Cs 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

WATCH: Pull the electronic energies from the output file.

Commands:

• grep "SCF Done:" 1.out will print to the screen every line in 1.out that contains the string “SCF Done:”

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

WATCH: Pull the frequencies from the output file.

Commands:

• grep "Frequencies" 1.out will print to the screen every line in 1.out that contains the string “Frequencies.”

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.

Geometry Absolute Energy (H) Binding Energy (kcal mol–1)
Monomer –76.4206270785
Dimer –152.854436576 –8.27

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. Below are binding energies of the water dimer as fully characterized at the listed methods and basis sets.

Ebind (kcal mol–1) NBasisa HF B3LYP B3LYP-D3BJ PBE0 PBEPBE $$\omega$$B97XD MP2
cc-pVDZ 48 –5.76 –8.27 –8.96 –8.26 –9.25 –8.03 –7.47
cc-pVTZ 116 –4.45 –6.09 –6.74 –6.27 –6.85 –6.18 –6.08
cc-pVQZ 230 –4.00 –5.32 –5.97 –5.63 –6.05 –5.49 –5.49
aug-cc-pVDZ 82 –3.91 –4.71 –5.35 –5.17 –5.28 –5.17 –5.26
aug-cc-pVTZ 184 –3.74 –4.57 –5.20 –4.98 –5.10 –5.01 –5.18
aug-cc-pVQZ 344 –3.72 –4.58 –5.21 –4.97 –5.11 –4.95 –5.09

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?