Quantumchemistry, NWI-MOL406, Computer exercises,
week 6 (21-Dec-2012) and week 7 (11-Jan_2013)
First a remark: next year we'll learn that DFT methods also use orbitals...
so it appears not that different from other methods (on casual inspection).
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 "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-multiparty 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) 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 AE 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 AE 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 AE = 0.62 eV). You'll observe good agreement, both
IP and AE, but there is something fishy here... what? how do you explain the variation of the AE
with basis set?
Week 3: crystal calculations
Carbon 3
NOTE: This doesn't work on the pc's in 00.026 because of some vague runtime error related to the Ubuntu release.
To do the exercises anyway: "ssh lilo". Please beware that you are not alone on lilo, so limit your memory to
1 Gbyte in the input file("memory 1000 mb"), and watch a bit the load of the systems (e.g. using "top"). The lilo has 8 cpu's
(or cores) and 8 Gbyte of memory in total. To leave lilo: type "exit".
Go the NWCHEM web site, go into the documentation, and fetch the example on the bandstructure of diamand
as calculated with "BAND". Run this example yourself (from the command line is easy, beware to set
appropriate paths for your files and scratch space). Make a little plot of the band structure (e.g. using xmgrace),
you can check the result on the NWCHEM web site. Why is the programme run twice?
Now adapt the input file, to calculate the band structure of a graphene layer (a 2-dimensional honeycomb lattice,
the lattice constant is 2.46 Angstroms). You have to mimic a real graphene by choosing a large distance
(say 7 Angstroms) between the layers. For the self-consistent part of your calculation, use an 8x8x1 mesh.
Make a plot going from M to Gamma to K (zoom in on the region around K). Also check your band structure
(you may find differences with "a" below in the empty states, but your calculation should be better).
Where is the "Fermi level"?
Links:
graphene a,
graphene b
Do not forget de-install before you log-off.
- Back to [ Quantum Chemistry ]
Last updated: 9-Jan-2013