Quantumchemistry, NWI-MOL406, Computer exercises,
We'll start with the "black font" exercises.
A severe test for DFT methods is the calculation of atomization energy (the energy
needed to break up a molecule into separate neutral atoms).
We'll calculate the atomization energy of H2, and next that of CH4.
Beware: use the online manual!
The hydrogen atom
As a pre-requisite, we first calculate the H atom.
As basis set take: 6-311G** (overkill, but we don't want to change it later).
Do an energy calculation for 3 DFT functionals and HF:
1. a plain LDA: Slater exchange and Perdew-Zunger (1981) correlation,
2. a good gradient corrected functional: PBE (Perdew, Burke and Ernzerhof),
3. a standard hybrid functional (DFT and HF mixed): B3LYP.
4. standard HF (you can do that from the "dft" mode, using the keyword "XC HFexch")
Note that this atom is an open shell system (look in the manual how to set multiplicity).
Make a list of the energies and a list with the eigenvalues... amazed?
The atomization energy of H2
Now we're going to calculate the atomization energy of H2.
Build the molecule, optimize the structure, and calculate the frequencies.
Do it for the 3 different DFT functionals above.
Make a list of total energies, bond lengths and frequencies.
The experimental bond length is 0.74 Angstrom,
The experimental atomization energy at 0 Kelvin is 103.5 kcal/mole.
Why did you have to calculate the frequencies?
Again reflect why you might be amazed...
The atomization energy of CH4
Now you have done this for H2, do it for a slightly more challenging molecule: CH4.
This example ought to show the superiority of 2 and 3 over 1.
Note: CH4 is a tetrahedral molecule.
Note: think about the spin-multiplicity of C.
Inspect the eigenvalues of the alpha and beta orbitals.
Note that in LSDA (local spin density functional) theory we have two densities,
one for the alpha and one for the beta spin.
The experimental atomization energy is 392.5 kcal/mole.
A chemical reaction
Now we'll do a simple chemical reaction. We'll use CH4, H2 and C2H4. So first repeat the exercise
of CH4 for C2H4 (this should be routine now). Again you'll need an experimental number (531.9 kcal/mole).
Next calculate the (T=0K) enthalpy change of: 2 CH4 → C2H4 + 2 H2, both experimental and for the
3 functionals. What do you observe? Do things get better for the LDA? Why?
The water dimer
If time is left, you can start on the water dimer. Use B3LYP and optimize.
Start from the model below [from Halkier et al., Theor. Chem. Acc. (1997) 97:150-157]
that is already good.
O -1.41084 0.063712 0
H -1.74764 -0.415688 -0.756952
H -1.74764 -0.415688 0.756952
O 1.48172 -0.066612 0
H 0.53235 0.055606 0
H 1.83806 0.82179 0
The two molecules are linked via a hydrogen bond. Its strength (at 0 K) is approximately
3.5 kcal/mole. Calculate the strength in two ways: just by subtracting the
energy (forget about zero-point contributions) of the fragments from that of the complex,
and by doing the same with BSSE (Basis Set Superposition Error, see note) correction
(see dft part of manual). What is happening here? How do energies change with this so-called
counterpoise correction method?
(you don't need to further optimize the isolated molecules, that
effect should be minute, but you can check).
Adding and removing an electron
Now we're goint to look a bit into ionization potentials and electron affinities.
First revist you results on H. In EXACT DFT the highest occupied KS orbital should give you
minus the IP. How large is the error here? Again you're confronted with a large self-interaction error.
A better way of calculating the IP is via `Slater's transition state', as in the werkcollege. However,
NWCHEM doesn't allow a non-integer number of electrons, so we can't explore that option. Instead we'll
do it the usual way, i.e. in a ``Delta-SCF'' fashion. So you do separate calculations on Li+, Li and Li-, and
obtain the IP and EA as total energy differences. You can use the `CHARGE' keyword to make ions.
Just use plain LDA.
Do a series of calculations for various basis sets:
3-21G 6-311G* 6-311G** 6-311+G*
Calculate both IP and EA as a function of basis set, and also keep track of the KS eigenvalue of
the highest occupied KS orbital. You'll note something specials for the negative ion.
Compare to the experimental numbers (IP = 5.39 eV, and EA = 0.62 eV). You'll observe good agreement, both
IP and EA, but there is something fishy here... what? how do you explain the variation of the EA
with basis set?
Do not forget de-install before you log-off.
- Back to [ Quantum Chemistry ]
Last updated: 14-Nov-2014