Lab 3
Details of a Calculation

Figs/Lab3/title.jpg

3.1  Overview

The major part of the first two labs in this course was spent learning how the graphical interface works and how to define the system to be calculated. Building molecules, submitting calculations and assigning molecular charge and spin multiplicity were discussed. Some aspects of visualization and analysis of the results were also considered.
Each of the electronic structure calculations performed in these first two labs have a number of parameters associated with them and several choices needed to be made before a calculation was run, choices that affected the quality of the results obtained. Until now, we have simply taken the default choices provided by the ADFINPUT program. In this lab we will investigate the two most important choices that must be made before a calculation can be run: the choice of density functional and the choice of a basis set. Before doing so, we will consider another choice that must be made, the choice of the core size. This choice is less important simply because it has less influence on the quality of the final results.

3.2  The Dissociation Energy

In this lab we will be comparing different ways to calculate molecular properties. In order to do that, we will need to have a property to calculate. A useful property for this purpose is the bond dissociation energy and the related properties the atomization energy and the heat of formation.

3.2.1  Definitions

The dissociation energy of a bond, De is the difference in energy between a parent molecule and the energy of the system once a particular bond is broken and the two pieces are separated by an infinite distance
Figs/Lab3/3_1.jpg
Figure 3.1: The dissociation of a C-H bond of benzene

De = (E(A)+E(B))- E(AB)
(3.1)
The atomization energy is the energy required to split a molecule into its constituent atoms
EA = Eatoms- Emolecule
(3.2)
The heat of formation of a molecule is the difference in energy between a molecule and the atoms that make up that molecule in their naturally occurring elemental form.
Hf = Eelements - Emol
(3.3)
Some or all of these definitions should be familiar to you. Note that we are neglecting several important contributions such as vibrational energies and temperature effects to keep things simple. The heat of formation is more complicated still because it can involve elements in the solid or liquid phase. We include it here for completeness.
All of these properties have in common that they are all essentially energy differences between one set of molecules and atoms and another set of molecules and atoms. We have discussed how to calculate the energy of an atom or molecule in section 1.5.1 in Lab 1. Thus, we are in a position to calculate dissociation energies, atomization energies and heats of formation.
Example
The dissociation energy of H2 = E(H2)-2E(H). The atomization energy of H2 is equal to its dissociation energy.
Example

De(C-H(benzene)) = E(benzene)-(E(phenyl radical)+E(H))
(3.4)

EA(benzene) = E(benzene)-(6E(C)+6E(H))
(3.5)
Example
The binding energy of an ethene molecule to a zirconocene catalyst is equal to the energy of the zirconocene-ethene complex minus the sum of the energies of the bare catalyst and isolated ethene
Exercise 3.1
Calculate the energy required to break the C-H bond in benzene and the atomization energy of benzene.
Note
Assume that all molecules and atoms are neutral in charge. The pieces formed in the dissociation/atomization of benzene are radicals. Be careful that you get the spin states correct and use an unrestricted calculation when you have unpaied electrons. Remember to clean up your structures with the happy face button after you build them.

3.3  Core Electrons

The idea of valence and core electrons is very important in chemistry. It allows the chemistry of the elements to be rationalized in a simpler way by assuming that only the highest energy valence electrons are chemically active while the lower energy core electrons are chemically inert. This idea underlies the periodic table. Elements with similar arrangements of valence electrons but different cores fall in the same column of the periodic table.
Computational chemists make use of the core-valence idea to simplify their calculations. The core electrons are less important for chemical problems and can either be
In a given calculation it is not always obvious which electrons should be put in the core and which should be called valence electrons. In a calculation electrons can be put into the core if doing so will make an insignificant difference to the calculated property of interest. The decision to put electrons in the core or valence will depend on the element in question, the property to be studied and the desired accuracy.
Example
The 1s electrons of C can almost always put put in the core unless very high accuracy is required.
Example
The 3p electrons of the first-row transition metals can usually be put in the core but sometimes should be in the valence, especially if higher accuracy is needed.

3.3.1  Core Electrons and ADF

In ADF the core electrons can be frozen. There are two ways to do this. We will consider the simpler but less flexible method here. On the right hand side of the ADFINPUT window amongst several options is a Core Type option (see figure 3.2) The drop-down menu has four choices: None, Small, Medium and Large. "None" means no electrons are frozen. The number of electrons frozen increases from the "Small" core to the "Medium" core to the "Large" core. Exactly what each of these means depends on the elements involved. For instance, for carbon, Small, Medium and Large all mean the same thing, a 1s core. For iron, Small means a 1s2s2p core while Medium or Large mean a 1s2s2p3s3p core. It is not possible to describe all the possible core choices in a few rules but in general lighter elements (2nd and 3rd row) have only one core choice (as well as no core) while the other elements have two choices, Small and Large. A few elements have all three choices.
Figs/Lab3/Core.jpg
Figure 3.2: The Core Type button.
The default choice is a Large core.
Exercise 3.2
Calculate the atomization energy of benzene with the "None" core option. Compare your result with the atomization energy of benzene that you calculated earlier.
Exercise 3.3
If you were interested in calculating the ionization potential of a 3d electron of iron, which core would you choose? What if you were removing a 3p electron? A 1s?

3.4  Basis Sets

When preparing to run a calculation, one of the first choices that must be made is the choice of a basis set. With ADF we will do calculations that use density functional theory (DFT). In a DFT calculation we attempt to calculate the orbitals and electron density of a molecule. Roughly speaking, a basis set is a set of mathematical functions available to describe the electron density and orbitals. In a calculation we determine what is the best possible way to combine these basis functions. By "combined" we usually mean added together. How do we find the "best possible way" to combine these basis functions, when we may have hundreds or thousands of them? That is one of the main computational tasks of any calculation and is one of the reasons so much computer time is needed.
There are many possible choices for basis sets. Usually more functions is better. A basis set can be thought of as the colours available to an artist who is painting a picture. If the artist has few colours (few basis functions or a small basis set) then she will not be able to paint a very accurate picture of the subject. If the artist has many colours (many basis functions or a large basis set) then she will have a better chance of painting an accurate picture of the subject. It is the same with basis functions. If we have many functions then we will have a better chance to correctly describe the orbitals and electron density of the molecule that we are interested in. Note that an intelligent choice of colours is just as important to the artist as the number of colours. An artist with millions of red, orange and yellow colours at hand will paint a less accurate picture of a tree than an artist who has a few browns and greens to choose from. The same idea holds for basis sets. Careful selection of basis functions can help enormously.
Most basis sets for molecular calculations are built of functions based on atomic orbitals centred on nuclei since these functions provide a good starting point. Each atom will have a set of functions at its location. The exact composition of that set will depend on a choice made by the user and on what element that particular atom is. Thus an atom of H might have a set of functions with s symmetry placed at its location while an atom of Mn should have a set of functions that include at least s, p and d functions associated with it. The combination of all these atomic basis sets is the overall basis set used in the molecular calculation. The number of basis functions used to describe the core electrons will be few or none depending on how the core electrons are treated.

3.4.1  Basis Sets and ADF

ADF uses atom-centred basis sets made up of Slater-type orbitals (STOs).
fSTO = Nrn-1e-zrYlm
(3.6)
The solutions of the Schrödinger equation for hydrogen are STOs so you may have seen the above equation before. In the equation, N is a normalization constant, n is the principle quantum number, z is a constant that is chosen somehow, r is the distance from the atom and Ylm is a spherical harmonic (angular function) that depends on the l and m quantum numbers. The Ylm determines whether a given basis function is s-type, p-type (x, y or z) and so on.
Figs/Lab3/Basis.jpg
Figure 3.3: The detailed basis set description in the output.
As was noted earlier, the basis set requirements of Mn are different to those of hydrogen and so each element has its own basis sets. In ADF, each element has several basis sets to choose from often including (but not limited to) SZ, DZ, DZP, TZP, TZ2P and QZ4P. Here S, D, T and Q stand for "single", "double", "triple" and "quadruple" and Z stands for "zeta" (z). Thus, SZ is called single-zeta and means one set of functions (optimized z from the above equation) describing each shell of valence electrons of the given atom. For H this means a single s function. For C this means a single s function and a set of three p functions (x,y and z). For Fe this means a set of d functions (x2-y2, z2, xy, xz and yz) for the 3d electrons and a single s function for the 4s electrons. DZ is double-zeta and means two sets of functions describing each shell of valence electrons of the given atom. Thus for H we have two s functions with different values of z, for C we have two s functions and two sets of p functions. TZ and QZ follow similarly. The P in TZP, TZ2P and QZ4P stands for "polarization." When an atom becomes part of a molecule the functions already mentioned may not be enough to describe the molecular orbitals well enough. Rather than adding more functions with different z values, it is often better to add functions with higher angular momentum (l quantum number) than those already present. These functions are the polarization functions. The benefit of these high l functions is their different angular behaviour (see http://www.orbitals.com/orb/orbtable.htm for pictures of many orbitals). The TZP basis for H therefore has three s functions and a set of p functions and the TZP basis for C has three s functions, three sets of p functions and a set of d functions. The TZ2P basis set has polarization functions of the next two higher l values. Thus TZ2P basis for H is the same as TZP but with an additional set of d functions and the TZ2P set for C is the same as the TZP set but with an additional set of f functions. The QZ4P basis set has two sets of polarization functions of the next two higher l values. Thus, the QZ4P basis set for H has four s functions, two p polarizations sets and two d sets. The QZ4P basis set for C has four s functions, four sets of p functions, two sets of d functions and two sets of f functions. The transition metal and f-element basis sets don't exactly follow these rules but are not far off.
If core electrons are frozen they are treated differently and are not described by basis functions. If the "None" option is chosen then the core electrons must also be described with basis functions. There is one function per shell in the SZ basis, two functions per shell in the DZ, DZP, TZP and TZ2P basis sets and three functions per shell in the QZ3P basis. It is also worth noting that frozen cores are not possible with the QZ3P basis at this time so the option "None" is always chosen with the QZ3P basis.
The default basis set used by ADF, and the choice we have used so far is DZ.
The detailed output of a calculation includes a lot of information about the basis set used in the calculation. To look at this information, open the OUTPUT program with SCM: Output. Select and output to look at. Once it is open, choose the Properties: Basis and Fit Functions option and you will find the basis set information (figure 3.3).
The basis sets of H and C are summarized in table 3.1. From this table we can evaluate that a calculation on H2 with a SZ basis set on each atom will have 2 functions, with a DZP basis set there with be a total of 10 functions and with a TZ2P basis set 22 functions will be used.
By using STOs, ADF is a little unusual. For reasons of computational economy, most computational chemistry packages make use of Gaussian type orbitals (GTOs) (with the form e-zr2).
Further information about basis sets can be found in several places on the web. Useful sites include http://en.wikipedia.org/wiki/Basis_set_(chemistry), http://www.ccl.net/cca/documents/basis-sets/basis.html and http://vergil.chemistry.gatech.edu/courses/chem6485/pdf/basis-sets.pdf. A very large database of GTO basis sets can be found at http://www.emsl.pnl.gov/forms/basisform.html.
H C
Basis Composition Tot. no. of functions Composition Tot. no. of functions
SZ 1s 1 1s1p 4
DZ 2s 2 2s2p 8
DZP 2s1p 5 2s2p1d 13
TZP 3s1p 6 3s3p1d 17
TZ2P 3s1p1d 11 3s3p1d1f 24
QZ4P 4s2p2d 20 4s4p2d2f 40
Table 3.1: Basis sets for H and C (valence part only).

3.4.2  Counting Basis Functions

The size of the basis set chosen will influence how long a calculation will take. It is therefore useful to be able to count the number of basis functions ahead of time to get an idea of the size of a calculation. For example, the number of functions in the H and C basis sets are listed in table 3.1. A molecule containing 10 C atoms and 22 H atoms that is described with a DZP basis set will require a total of 10 ×13 + 22 ×5 = 240 basis functions.
Exercise 3.4
While it is possible to chose different types of basis set for different atoms in a molecule we will use the same type for all atoms in a molecule for the moment. Make up a table listing the six types of basis set and giving how many functions will be used in a calculation on the N2 molecule and a calculation of a molecule of benzene.

3.4.3  Calculations with Different Basis Sets

Figs/Lab3/3_2.jpg
Figure 3.4: The basis set button.
The simplest method to choose which basis set you are using is through the same "main options" window that was used to choose the core size. Just below the button used to make this choice should be a button with DZ on it. Pressing the control to the right of the DZ should give a menu of basis set choices. This choice will apply to all atoms in a molecule.
Exercise 3.5
Do this exercise in groups of three. Between the three of you, calculate the atomization energy of benzene with the six different basis sets considered so far. Use the default core size in each case. Make up a table listing the following details of these calculations for each basis set: Number of basis functions, calculated energy of benzene, atomization energy of benzene and calculation time. Plot atomization energy and calculation time vs number of basis functions.
Exercise 3.6
What do you expect will happen to the three properties graphed in the previous two exercises if the basis set size is increased further?

3.5  Choice of Functional

In this course we have chosen to utilize density functional theory as our quantum mechanical method. However, once DFT has been chosen, the approximations in the theory still have not yet been specified completely. In DFT the energy of a system is written as a functional of its density.
E[r] = T[r] + Eext + EC[r] + EXC[r]
(3.7)
where r is the electron density, T is the kinetic energy, Eext is the energy due to an external potential such as the charge on any nuclei present, EC is the Coulomb interaction between the electrons and EXC is the so-called exchange-correlation contribution that includes everything that is missing in the previous two pieces. The exact form of EXC is not known. Many attempts have been made to obtain EXC and each of these attempts define a different choice of functional and a different choice of approximation. The question is, which functional should be chosen for the problem at hand?
Many, many many possible functionals have been proposed, tested and then published in the literature. Each is justified somehow and has its own strengths and weaknesses. So far, none of the available functionals is universal in the sense that it is very good for all possible applications. Thus, at present generally it is necessary to decide what is to be calculated and then choose the most appropriate functional of those available.
Functionals come in six rough classes
Usually corrections are added to an LDA to give a GGA or meta-GGA. The exact exchange of a hybrid functional is usually added on top of a GGA and LDA. GGAs often come in two parts a "correlation" and an "exchange" contribution that is related to but different from the exact exchange of a hybrid correction.

3.5.1  Functionals and ADF

Figs/Lab3/3_3.jpg
Figure 3.5: Choosing a functional.
In ADF you have the choice of over 50 DFT functionals including examples from all of the classes just mentioned except OEP. Up until this point we have been using the default functional which is a type of LDA.
The functional to be used in a given calculation can be chosen by selecting the XC Functional option from the Model menu. To choose a functional of the LDA, GGA, Hybrid or "Other" type click and the drop-down menu next to "xc functional in SCF". The meta-GGAs are accessed by clicking any of the two buttons. Making use of the menu will choose only one functional. Selecting one of the meta-GGA buttons will specify that all available meta-GGAs will be calculated (figure 3.5).
Note that no electrons can be frozen when applying a Hybrid functional. So, if we are running a calculation with a hybrid functional then we must choose a "NONE" core.
The total energy obtained from a GGA, hybrid or "other" calculation can be found in a similar place to the LDA energy (see section 1.5.1 in Lab 1). Often the X (exchange) and C (correlation) parts of a GGA are separate and their contributions are listed individually (figure 3.6).
Figs/Lab3/3_4.jpg
Figure 3.6: The energy from a GGA calculation.
The meta-GGA energies are obtained from a different place. In order to find them, the detailed output must be opened through the SCM: Output option. Once the OUTPUT program has started go to the energies section through the menu item Properties: Bonding Energy Decomposition. Scroll down a little further and you will come to a long list of energies including all available meta-GGAs (figure 3.7)
Figs/Lab3/3_5.jpg
Figure 3.7: The energies from a meta-GGA calculation.
Exercise 3.7
Work in groups of three for this exercise. Between the three of you, calculate the atomization energy of benzene with six different functionals (except LDA). At least two of the functionals should be GGAs, at least one a hybrid functional and at least one a meta-GGA. Use a DZP basis set. For each functional, tabulate the atomization energy, and calculation time. Include your LDA/DZP data from the previous section in your table. Graph you results vs choice of functional. Order the functionals so that LDA is first, followed by the GGAs, meta-GGAs, hybrids and others (if any) in that order.
Exercise 3.8
What trend in the calculated properties do you observe as the choice of functional goes from LDA to GGA to meta-GGA to hybrid to other?

3.6  What Works Best?

From a combination of experimental and theoretical data, it can be discerned that the atomization energy of benzene is about 5501 kJ/mol.
Exercise 3.9
Work in pairs for this exercise. Each pair will choose a functional (not LDA) and calculate the atomization energy of benzene with that functional using the TZP and QZ4P basis sets. Evaluate the error in you results (E(yours) - 5501 kJ/mol) and report back to the class. Make a table of the results and find the best functional for calculating the atomization energy of benzene of those chosen. Is it necessary to use a QZ3P basis set to get an accuracy of better than 20 kJ/mol with the best functional or is TZP enough? Note that 1 hartree of energy is equivalent to 4625.5 kJ/mol.

3.7  Projects

If you are a graduate student, do project 3.7.1. If you are an undergraduate, do project 3.7.2, unless you really want to do project 3.7.1. Before you decide on exactly what you will do for the second project, consult your T.A. for advice. No two people can look at the same functional or journal article.

3.7.1  Definitions of Functionals

As has already been noted, functionals can be defined in many ways. Choose one of the functionals available in ADF. Find the original reference of that functional. Study that reference and then briefly describe how the researchers derived that functional. Your description should be short, only one or two paragraphs. Don't go into too much detail but just describe in a general way how the functional was defined.
You can choose from:
LDA
Describe the VWN parameterization.
GGA
A GGA usually consists of an exchange and a correlation correction to the LDA. Choose an exchange correction or a correlation correction. For example the BLYP functional is made up of a B piece and a LYP piece. If you are not sure consult with your T.A.
Other
Describe one of the model functionals.
Don't choose a hybrid or meta-GGA functional.
The references for the LDA, GGAs and Other (Model) functionals can be found in the ADF documentation at http://www.scm.com/Doc/Doc2006.01/ADF/ADFUsersGuide/page82.html.

3.7.2  Tests of Basis Sets and Functionals in the Literature

When a functional is proposed it is generally not clear how good that functional will be. The accuracy of a given functional is usually demonstrated through example calculations where theoretical results are compared with experimental data. The field of chemistry includes an enormous number of types of molecules, bonding situations and properties that might be interesting to calculate. The task of testing each functional against the possible applications is an enormous one. Many papers have been published that have tackled a small part of this problem.
In this project you will look at one of them. Choose one paper from the list below. All of these papers describe calculations where different basis sets and/or functionals have been tested somehow. Read the relevent parts of the paper and in a couple of paragraphs outline
Note that some of the papers listed below contain other material that does not deal with the tests. You do not need need to say anything about these sections. We are only interested in the parts where basis sets and/or functionals were tested.
Here is the list of papers. Electronic versions of all of these papers can be accessed through the University of Calgary website.
  1. Erik van Lenthe and E. Jan Baerends, J. Comp. Chem. 24, 1142-1156 (2003). This paper describes how some of the basis sets used by the ADF program were created and also tests those basis sets. If you want to look at this paper you can describe one of two things. Either how the exponents of the basis functions were optimized or the tests looking at the quality of the basis sets. If you choose the first option you will find that most of the details are in fact in the first few pages of a different paper (J. G. Snijders P. Vernooijs and E. J. Baerends At. Data Nucl. Data Tables 26 483-509 (1981).)
  2. Nick X. Wang and Angela K. Wilson J. Phys. Chem. A. 107, 6720-6724 (2003).
  3. Jianmin Tao and John P. Perdew J. Chem. Phys. 122, 114102 (2005). Do not worry about section III.
  4. Marcel Swart, André R. Groenhof, Andreas W. Ehlers and Koop Lammertsma J. Phys. Chem. A. 108, 5479-5483 (2004).
  5. Nathan E. Schultz, Yan Zhao and Donald Truhlar J. Phys. Chem. A. 109 4388-4403 (2005). This paper contains lots of results. When you are considering which functionals the authors judged to be the best, concentrate on the Conclusions section.
  6. Miriam M. Quintal, Amir Karton, Mark A. Iron, A. Daniel Boese and Jan M. L. Martin J. Phys Chem A. 110 709-716 (2006). These people compare their DFT results with other calculations. Don't worry about what the other calculations (labelled CCSD(T)) mean. You can assume that the CCSD(T) numbers are very accurate. Also do not worry too much about any discussion of the "% HF exchange".
  7. Sergey N. Maximoff, Juan E. Peralta, Verónica Barone and Gustavo E. Scuseria J. Chem. Theory Comput. 1, 541-545 (2005).
  8. Filipp Furche and John P. Perdew J. Chem. Phys 124, 044103 (2006). This paper also includes lots of results. When you are considering which functionals the authors judged to be the best, concentrate on the Conclusions section. Don't worry about the theory section (II) or the appendix.
  9. Charles W. Bauschlicher Jr. Chem. Phys. Lett. 246, 40-44 (1995).


File translated from TEX by TTH, version 3.85.
On 15 Sep 2009, 07:08.