SC200 - Report 3 - Ragnhild Blikberg and Thierry Matthey

UiB : MatNat : Informatikk : Undervisning : Kurs : SC200 : SC200 Thierry Matthey

SC200 - Report 3 - Ragnhild Blikberg, Thierry Matthey

Computational Chemistry

[Problem Description] [Introduction] [The Model] [Num2D] [Gaussian] [Conclusion] [Constants and Units] [References]

Introduction:

The purpose of this exercise is to use some packages from computational chemistry. The first package, Num2D, is quite simple and only applicable to diatomic molecules. We will also look at the much more powerful package, Gaussian, which is used worldwide in chemical physics.

Its important to understand the interaction between molecules and electromagnetic radiation. What happens with a molecule which absorbs fotons depends on the wavenumber:

Microwaves Rotation

Infrared Vibration

Visible light Electronic excitation Colors

UV Excitation, ionization Ozonlayer

IR spectroscopy is an important technique in chemical analysis of an unknown probe. The reason is that different functional groups absorb infrared light only in given wavenumber intervals. IR spectroscopy is also much used to get information about the bonding in a molecule. The stronger a bonding is, the shorter bond lengths and higher frequencies of vibration.

All constants used in the report are given in [Constants and Units].

The Model

[Problem Description] [Introduction] [The Model] [Num2D] [Gaussian] [Conclusion] [Constants and Units] [References]

The molecular Schrödingers equation, apart from electron spin effects, is

where

is the Hamiltonian and E is the total energy. The different terms in the Hamiltonian are given by

The molecule consists of N_e electrons and N_n nuclei. The electronic position coordinates are designated r₁,r₂ ..., while the nuclei position coordinates are designated R₁,R₂ ... The nuclei masses are given by M₁,M₂ ... and m_e is the electronic mass. Z₁,Z₂ ... are the charges of the nuclei.

The wavefunction that satisfies the Schrödinger equation must be a function of both the electron position coordinates and the nuclear position coordinates, and this PDE is not separable. This problem will be dealt with through the Born-Oppenheimer approximation.

First we will restrict ourselves to diatomic molecules, i.e. N_n=2. Then we need only one coordinate R to specify the separation of the nuclei. As mentioned, we will use the Born-Oppenheimer approximation. This leads to a separation of the molecular Schrödinger equation into one part for the electronic wavefunction and one part for the nuclear motions, which is the Schrödinger equation for vibration and rotation. Thus, we approximate the molecular wavefunction by a product

The essential element in the approximation is that applying the nuclear kinetic energy operator to the electronic wavefunction yields zero. The physical idea is that the light, fast moving electrons readjust to nuclear displacements instantaneously. This gives the following electronic Schrödinger equation where R is now a parameter

where

The nuclei part of the wavefunction must then satisfy the following equation

The Hamiltonian consists of a kinetic energy operator for the nuclear position coordinates, the repulsion potential between the nuclei, and an effective potential, in the form of E_e, that gives the energy of the electronic wavefunction as it depends on the nuclear position coordinates. For a diatomic molecule this is a one-dimensional equation. For molecules with more than two nuclei, the dimension of this equation will increase (see below). The combination of E_e and V_n-n, yields the effective potential for molecular vibration. It is the potential energy for the nuclei in the field created by the electrons.

Solution of the electronic Schrödinger equation:

Density functional theory
Expansion of orbitals in known basis

Solution of the Schrödinger equation for the nuclei motion:

Harmonic approximation of E_e+V_n-n
(3N_n -6) dimensional problem

Num2D

[Problem Description] [Introduction] [The Model] [Num2D] [Gaussian] [Conclusion] [Constants and Units] [References]

First, nuclear motion in diatomic molecules will be studied neglecting interactions with other molecules. Num2D is a numerical 2D package for solving the Hartree-fock and Dirac equations for diatomic molecules. We downloaded the package at http://www.csc.fi/~laaksone/num2d.html.

Problem 1

We chose the diatomic molecules H₂⁺ and H₂. Using the Num2D code, we calculated the effective potential energy for molecular vibration U(R) = E_e(R)+V_n-n(R) from the molecular Schrödinger equation for different values of R (nuclear separation) with dR = 0.08. The results are given in h2.txt and h2p.txt.

Problem 2

From [2] we found that the potential energy approximated by Taylor's formula is

If and only if R₀ is a point of stable equilibrium then U'(R₀) = 0, and the function gets a parabola i.e. the harmonic approximation

where U(R₀) is the potential energy of equilibrium and K = U''(R₀) is the curvature of the parabola. To approximate K, one can put the results from Num2D into the central difference formula

Figure 1: H₂ molecule.

For H₂ (Figure 1), R is plotted (Figure 2, blue line) against the effective potential U(R) of problem 1. We have chosen dR = 0.08 and found R₀ = 1.36 a₀, U(R₀ - dR) = -1.131 E_h and U(R₀ + dR) = -1.133 E_h. From these results it was possible to approximate U''(R₀) and we got

Figure 2 also contains approximation by the harmonic potential given above. This was plotted by the MatLab program plot_h2_harm.m.

Figure 2: U(R) vs. R for H₂.

Figure 3: H₂⁺ molecule.

For H₂⁺ (Figure 3), R is plotted (Figure 4, blue line) against the effective potential U(R) of problem 1. We have chosen dR = 0.08 and found R₀ = 2.00 a₀, U(R₀ - dR) = -0.602307 E_h and U(R₀ + dR) = -0.602305 E_h. From these results it was possible to approximate U''(R₀) and we got

Figure 4 also contains the approximation by the harmonic potential given above and the analytic (exact) solution (green line) from [1]. This was done by plot_h2p_harm.m.

Figure 4: U(R) vs. R for H₂⁺.

For Figure 2 and 4 the harmonic approximation has a high accuracy close to the equilibrium state (R₀) while for R >> R₀ the approximation fails entirely. From Figure 4 we can conclude that the numerical solution is very close to the analytical (exact) solution.

Problem 3

We will now calculate the different energy levels the system can achieve by the nuclear Schrödinger equation:

The kinetic energy operator is

where M₁ and M₂ are the masses of the nucleis, and

where mu is the reduced mass. Inserting the kinetic operator (T_n) and the harmonic potential of the effective potential

into the Schrödinger equation yields

Reorganizing yields

The energy levels corresponding to the harmonic potential according to [2] is given by

where n is the vibration quantum number and the vibration frequency is given by

Table 1 shows the calculation of omega for H₂ and H₂⁺.

Table 1: K, reduced mass and frequency for H₂ and H₂⁺

If a closer look at the calculations is needed, look here.

This means that the allowed energies of the system is not continues, but quantized. The energy levels corresponding to the harmonic potential is plotted with plot_h2_harm_elevels.m and plot_h2p_harm_elevels.m and given in Figure 5 a) H₂ and b) H₂⁺.

Figure 5 a) and b): Energy levels corresponding to the harmonic potential for H₂ (left) and H₂⁺ (right)

From the Figure 5 a) and b) we observe that the energy levels corresponding to the harmonic approximation are equally spaced.

One can approximate the potential by functions with higher accuracy than the harmonic potential. We chose the basis functions {x², x, 1, x^-2, e^x, e^-x} and solved a linear system with 6 equations and 6 energy values obtained by Num2D. For H₂ the solution of the linear system gave us the function

and for H₂⁺ the function

These functions will give the effective potential in [eV].

A package SLEIGN2 at http://www.math.niu.edu/~zettl/SL2/ can be used for solving the Sturm-Liouville system (nuclear Schrödinger equation)

By this package we calculated the energy levels using our new approximations of the effective potential. The energy levels are plotted with plot_h2_elevels.m and plot_h2p_elevels.m in Figure 6 and 7.

Figure 6: Energy levels for H₂

Figure 7: Energy levels for H₂⁺

With the more accurate approximation of the potential energy, we get an asymptotic close distribution of the energy levels when R -> infinity.

We remark that for low energy levels, the harmonic and the new approximation correspond. For higher energy levels (large R) they differ.

Gaussian

[Problem Description] [Introduction] [The Model] [Num2D] [Gaussian] [Conclusion] [Constants and Units] [References]

Problem 1

In organic chemistry, ethane, ethene and ethyne are the classic examples of molecules with simple, double and triple carbon-carbon bonding. Benzene is the classic example of an aromatic molecule where it is impossible to assign an integer bond order. Figure 8, 9, 10 and 11 show the molecule structure and bonding of ethane, ethene, ethyne and benzene respectively.


Figure 8: Ethane, C₂H₆, simple CC bonding	Figure 9: Ethene, C₂H₄, double CC bonding

Figure 10: Ethyne, C₂H₂, triple CC bonding	Figure 11: Benzene, C₆H₆, 1-2 bonding

For the visualization of the molecules we used the raytracer from BOOGA with BSDL as scene description language. The BSDL files include a file to define common properties of the molecules, i.e. color, radius etc.. Additionally we generated VRML files of each molecule, to take a closer look at them directly via a web-browser (If your web-browser does not support VRML, we suggest using VRwave. At UiB, this is installed at /net/hstud/sparc/i291/vrwave-0.9/vrwave).

Problem 2

There are different ways of assigning intermediate bond orders. We want to do this by using the CC-stretch frequency as an index. Gaussian is a chemistry package which is widely used to solve quantum chemistry problems numerically. We have used Gaussian 98 to calculate the CC-stretch frequency for ethane, ethene and ethyne. For ethane, we started the program by the command

g98.job < etan.inp

This resulted in an output file etan.out. Similarly for the other molecules. The input-, output-, BSDL- and VRML-files are given in Table 2 below.

Molecule	Input	Output	BSDL-file	VRML-file
Ethane	etan.inp	etan.out	etan.bsdl3	etan.vrml
Ethene	eten.inp	eten.out	eten.bsdl3	eten.vrml
Ethyne	etyn.inp	etyn.out	etyn.bsdl3	etyn.vrml
Benzene	benzen.inp	benzen.out	benzen.bsdl3	benzen.vrml

Table 2: Input-, output-, BSDL- and VRML-files for ethane, ethene, ethyne and benzene.

With the visualization tool Molden we measured and found the different CC-distances and CC-stretch frequencies for each molecule. For ethane the command for using Molden was

molden35 etan.out

The results for all the molecules are given in Table 3.

Molecule	CC-stretch frequency [cm^-1]	CC-distance [Angstrom]	Bond order
Ethane	1004.343994	1.5351	1
Ethene	1522.09041	1.339	2
Ethyne	2233.886963	1.20	3
Benzene	1234.67004	1.38617	1.4893^*

Table 3: CC-stretch frequency in cm^-1 and Angstrom, and bond order for ethane, ethene, ethyne and benzene. ^* interpolated value from problem 3.

From Table 3 we can observe that the higher (shorter) frequency (distance) for the molecules (between the C atoms), the higher bond order.

Problem 3 and 4

From the results in Table 3 we interpolated with a quadratic polynomial through the points bond order vs frequency for ethane, ethene and ethyne. With the interpolant we computed the intermediate bond order of benzene, interpolate.m. The curve is plotted in Figure 12. We also did the interpolation for the CC-distance (bond length) for all molecules, Figure 13.

Our approximated bond-oder value for benzene is 1.4893, which is quite close to the theoretical value 1.5. For the CC-distance (bond length) we did not get such a good approximation (1.7292) as for the CC-stretch frequency.

Figure 12: CC-frequency vs bond-order

Figure 13: CC-distance vs bond-order

Conclusion:

[Problem Description] [Introduction] [The Model] [Num2D] [Gaussian] [Conclusion] [Constants and Units] [References]

In the first part of this exercise, we studied the nuclear motion in the diatomic molecules H₂ and H₂⁺. Using the numerical 2D package Num2D, we calculated the effective potential energy for molecular vibration from the molecular Schrödinger equation for different values of nuclear separation. As expected the distance (bond lengt) for the ground state for H₂⁺ is larger than for H₂. We then approximated the potential energy by a harmonic approximation. We conclude that the harmonic approximation has a high accuracy close to the equilibrium state, while for nuclear distances much larger than the equilibrium state, the approximation fails entirely. For H₂⁺ we compared the numerical solution from Num2D with the analytical solution [1] and conclude that numerical solution is very close to the exact solution.

From the Schrödinger equation it was possible to compute the energy levels of the system using different representations of the potential energy. For the harmonic approximation, we calculated this analytically, and got uniformly spaced energy levels. Using a more accurate approximation and solving the Sturm-Liouville system by the package SLEIGN2, the energy levels became non-uniformly spaced. The first levels corresponded with the levels achieved by the harmonic approximation, while the levels quickly became closer and closer, approaching the potential energy when the nuclear separation is large i.e. -> infinity.

In the second part of the exercise we studied the hydrocarbons ethane, ethene, ethyne and benzene. First we generated pictures of each hydrocarbon illustrating the bonding. With the help of the packages Gaussian and Molden we calculated the CC-stretch frequency and CC-distance for all molecules. By interpolation of the bond orders of ethane, ethene and ethyne we calculated the bond order for benzene. Using the CC-stretch frequency we got 1.4893, close to the theoretical value 1.5, but with the CC-distances we got a less accurate value 1.7292.

We conclude that the higher (shorter) frequency (distance) for the molecules (between the C atoms), the higher bond order. To compute the bond order the CC-stretch frequency gives better accuracy compared to CC-distance interpolation.

Constants and Units:

[Problem Description] [Introduction] [The Model] [Num2D] [Gaussian] [Conclusion] [Constants and Units] [References]

Constants and Units	Symbol	Formula	Value
Bohr radius	a₀		5.29E-11 m
Ångstrom	Å		10E-10 m
Electron rest mass	m_e		0.109E-31 kg
Atomic mass unit	u		1.6606E-27 kg
Electron charge	e		1.602E-19 C
Electron Volt	eV		1.602E-19 J
Planck's constant	h		6.626E-34 Js = 4.136E-15 eVs
	h	h/(2pi)	1.055E-34 Js = 6.582E-16 eVs
Hartree energy	E_h	h²/m_ea₀²	4.36E-18 J = 27.211 eV
Frequency	1 [cm^-1]	corresponds to	29 979 245 800 [Hz]

References:

[Problem Description] [Introduction] [The Model] [Num2D] [Gaussian] [Conclusion] [Constants and Units] [References]

[1] 'Physics of Atoms and Molecules', B. H. Brandsen and C. J. Joachim, Longman, 1983.

[2] 'Modern Physics', Serway, Moses, Moyer, 1989.

Last update by matthey@ii.uib.no and ragnhild@ii.uib.no