Chem350: Thermochemistry using Gaussian. In this write-up I will quickly go through the features needed to obtain thermochemical data using gaussview/gaussian09. I will not give you very detailed instructions. You will have to find your way using the graphical interface, which is quite straightforward to use. I encourage you to draw some molecules of your own, and see if you can manipulate the various tools in Gaussview. Below I will give you some instructive examples that illustrate many of the capabilities of Gaussview/Gaussian09. In lab don't hesitate to ask questions. The primary purpose of this section of Chemistry 350 is to show how the formalism we developed in class directly translates into simple formulas that involve quantities that can be extracted from routine quantum chemistry calculations: Optimized geometries and moments of inertia, molecular mass, vibrational frequencies and electronic energies. Together with this lab I created a Matlab file that requires these quantities as inputs and it then proceeds to calculate thermodynamical properties. We will check that the results are exactly the same as in the Gaussian program. Moreover, this gives you explicit access to the proper formulas and unit conversions. A. Calculations on a simple example: single point, geometry optimization and vibrational frequencies of H2O. Let us run through a very simple example first. You can create the water molecule using the builder (pull down the menu), selecting Oxygen from the ‘atoms’ menu. Let us first do a calculation on the guessed structure. Measure the bond angles and distance first (using inquire, ?). Let us do a Hartree-Fock / 6-31G(d,p) calculation. For this click Calculate and select Jobtype=Energy, Method=Hartree Fock, Basis is 6-31G, and on the right hand options select from the between brackets options (d,p), by pulling down the appropriate menus. This specifies the atomic orbital (AO) basis set. Then Submit the calculation (see lower-left button). You will have to save your file, for example as h2o. Follow the directions given by the program. In a little while you get back the result, and will be asked if you want to see it. Look at results (open the file h2o.log) and look under summary. Here the main specifics of the calculation are summarized and you can find the total energy. Next we will optimize the geometry, selecting jobtype=Optimization. Measure the geometries again. You should find the results as listed in the Table below. Table I. Results for water Property / units?

Guess geometry

R(OH) A(HOH) HF-6-31G(d,p) Energy

0.96 109.471 -76.02255367

1

HF – 6-31G(d,p) optimized geometry 0.94319 105.921 -76.02361493

At the optimized geometry we can now calculate vibrational frequencies. Start from the geometry in the h2o output file and select frequency under calculate->Jobtype. After the result returns you can view the vibrational frequencies (select under results). You can view normal modes and infrared and raman spectra. The frequencies for the normal modes in water are as follows. Normal Mode Bending Symmetric stretch Asymmetric stretch

Frequency (cm-1) 1770(.5) 4145(.87) 4262(.48)

This same calculation can be done in one shot. Start anew by rebuilding H2O from scratch, and select Opt+Freq under the Jobtype menu. It first does the optimization, and then calculates the vibrational frequencies. Vibrational frequencies are only meaningful at the optimized geometry! You can also look in the output file. The vibrational frequencies are all listed near the end of the file. This section also lists their symmetry. All of the vibrational motions have a definite symmetry with respect to planes of symmetry and rotational axes in the molecule. Normal modes of A1 symmetry preserve the symmetry of the molecule. Normal modes of other symmetry-type break the symmetry of the molecule. You can investigate by visualizing normal modes, under summary. B. Calculations on a slightly more advanced example: trans-butadiene. Some more options and features in Gaussian. I will run you through one more example doing similar things. Part of the information is similar as for water. Use the gv builder to construct trans butadiene. You only need to build the carbon chain using the appropriate type of carbon atom. The hydrogens are added automatically. Go to View to examine different ways to see the structure. In particular switch on Symbols to see all atom types. Rotate and translate the molecule, using the mouse. If you press the clean button (broom) in the builder window, you will get a first rough guess for the structure. It is also useful to click the symmetrize button. It can be useful to use the pointgroup option under edit. Your molecule should have C2h symmetry. Save your input file, e.g. t_c4h6. H H

c3 c1

H

H c4 H

c2 H

2

Next optimize the structure and calculate frequencies at the HF/3-21G level of theory. To this end open the calculate menu and specify options as before. - Measuring geometrical features. Click on the inquire button in the builder menu (question mark). This will give you geometrical information. By clicking two atoms you get the bond distance. By selecting three atoms you find the angle made by these atoms. Finally by clicking 4 atoms in sequence you get the dihedral angle (in "4321" format). This is the torsional angle around atoms 2 and 3. The dihedral angle is measured rotating the fourth atom onto the first atom along the 3-2 bond axis. Clockwise angles are positive. Now find the C1-C2, C2-C3 and C3-C4 bond lengths. You clearly see the double bonds in the geometrical structure. However, also C2-C3 is substantially shorter than typical C-C single bonds (quickly construct and optimize ethane if you wish to confirm this). Compare terminal and center CH bond lengths. Also measure various bond angles and comment on your findings. What is the C1C2C3C4 torsional angle around the C2-C3 bond? Now measure it by clicking the appropriate atoms. It is of interest to look at the geometry in the original guessed structure, before optimization. Here you see there is no information on conjugation. Instead the program guesses standard single and double bonds. Under the Results menu you will find a number of options that allow you to analyse the output of your calculation in an easy way. - Changing geometries. Let us transform trans-butadiene into cis-butadiene. To this end click on dihedral in the builder menu, and select four atoms that specify the desired dihedral angle (in '4321' format). This opens the dihedral smartslide. Change the dihedral from 180˚ to 0˚. Save this new structure as cis-c4h6.com and perform an optimization of the geometry. Now you can analyse the structure (and energetic and orbitals) of cis-butadiene. You can similarly change angles and bond distances by clicking the respective buttons in Gaussview. What is the symmetry of Cis-butadiene? - Vibrational frequencies and viewing normal modes. It is very easy to obtain vibrational frequencies in gaussview. Under the jobtype menu specify OPT+FREQ. Infrared intensities are calculated always, Raman intensities sometimes by default. Calculate the vibrational spectrum of trans-butadiene. Upon completion of the calculation you can display by going to Results/Vibrations. This shows a list of the calculated vibrational frequencies. By clicking any one of the frequencies and clicking start you can view an animation of the vibrations. Optionally you can also display the displacement vectors that characterize the normal mode. Set the Frames/cycle and Displacement scales to your preference. By clicking spectrum you see the complete vibrational spectrum displayed. Characterize the various vibrations in transbutadiene. Where are the CH stretches and how many are there? What is the difference

3

between them? Where are the various C-C stretches? Can you clearly distinguish C1-C2 stretch and C2-C3 stretches? Characterize the various bending modes, or rotations around the various bonds. - Charges. For any calculation you can view the charges in the molecule. Trans-butadiene is not particularly interesting of course. You might create an alcohol or a carbonyl compound to see something more interesting. Play around with the options to obtain a rendering of the charges that you prefer. The above gives you an overview of some of the most widely applicable options in Gaussview/Gaussian09. C. Investigating Thermochemistry and checking Statmech formulas in Gaussian. By optimizing a geometry and calculating vibrational frequencies one obtains all of the information needed to calculate the partition function of a molecule in the gas phase (ideal gas / rigid rotor / harmonic oscillator approximation). In order to achieve accurate results that can be compared to experimental values one needs to do accurate electronic structure calculations. For many practical purposes a Density Functional Calculation using the B3LYP functional and a cc-PVTZ basis set is sufficiently accurate. For a molecule that has lone pairs or internal hydrogen bonds a bigger basis set may be required. We will first do calculations on a linear molecule (CO2) and on the water molecule. I will show you where you can locate the thermochemistry output in the Gaussian output file. Moreover , we can run a little Matlab program to do the statmech part of the calculations ourselves. This way we can check what is done in Gaussian, and if things agree we have a compact representation of all the required formulas and all of the unit conversions through the Matlab file. i) Thermochemistry of water (non-linear molecule) Create the water molecule. Under jobtype select OPT+FREQ, while under Method select a ground state DFT calculation. Select the B3LYP functional (default) and the cc-pVTZ basis set under the appropriate menu. Then save the file and submit the calculation. When  you  go  under  Results/Summary  you  will  find  the  total  electronic   energy:     E(RB3LYP) -76.45983962 (in Hartree) This  is  one  piece  of  information  needed  for  thermochemistry.    If  we  then  click  view file,  the  output  file  is  opened.  If  you  scroll  the  file,  near  the  end  you  find  the   information  on  the  vibrational  frequencies.  For  water  it  looks  like  this:    

4

Harmonic frequencies (cm**-1), IR intensities (KM/Mole), Raman scattering activities (A**4/AMU), depolarization ratios for plane and unpolarized incident light, reduced masses (AMU), force constants (mDyne/A), and normal coordinates: 1 2 3 A1 A1 B2 Frequencies -- 1638.7752 3803.1701 3903.6257 Red. masses -- 1.0831 1.0447 1.0819 Frc consts -- 1.7139 8.9032 9.7133 IR Inten -- 69.5441 3.2440 40.9209 Atom AN X Y Z X Y Z X Y Z 1 8 0.00 0.00 0.07 0.00 0.00 0.05 0.00 0.07 0.00 2 1 0.00 -0.43 -0.56 0.00 0.59 -0.39 0.00 -0.56 0.43 3 1 0.00 0.43 -0.56 0.00 -0.59 -0.39 0.00 -0.56 -0.43 For  us  the  important  part  concerns  the  Frequencies,  which  has  the  harmonic   vibrational  frequencies  in  wave  numbers.       Frequencies -- 1638.7752 3803.1701 3903.6257 Just  below  this  section  on  frequencies  in  the  output  file  we  find  all  the  information   on  thermochemical  quantities.  This  section  of  the  output  file  looks  like  this      -------------------   - Thermochemistry ------------------Temperature 298.150 Kelvin. Pressure 1.00000 Atm. Atom 1 has atomic number 8 and mass 15.99491 Atom 2 has atomic number 1 and mass 1.00783 Atom 3 has atomic number 1 and mass 1.00783 Molecular mass: 18.01056 amu. Principal axes and moments of inertia in atomic units: 1 2 3 Eigenvalues -- 2.21189 4.15886 6.37075 X 0.00000 0.00000 1.00000 Y 1.00000 0.00000 0.00000 Z 0.00000 1.00000 0.00000 This molecule is an asymmetric top. Rotational symmetry number 2. Rotational temperatures (Kelvin) 39.15821 20.82634 13.59554 Rotational constants (GHZ): 815.92565 433.95112 283.28541 Zero-point vibrational energy 55898.9 (Joules/Mol) 13.36016 (Kcal/Mol) Vibrational temperatures: 2357.83 5471.91 5616.44 (Kelvin)

5

In sequence we find information on the atomic and molecular masses, the moments of inertia tensor (both the eigenvalues and the eigenvectors), rotational temperatures, and the symmetry number. Finally, vibrational temperatures for each normal mode. This is all the information needed to process the data and calculate heat capacities, entropy, Internal energy, Enthalpy, Gibbs free energy. The Gaussian program does it for you of course, and the results are listed as follows: Zero-point correction= 0.021291 (Hartree/Particle) Thermal correction to Energy= 0.024126 Thermal correction to Enthalpy= 0.025070 Thermal correction to Gibbs Free Energy= 0.003649 Sum of electronic and zero-point Energies= -76.438549 Sum of electronic and thermal Energies= -76.435714 Sum of electronic and thermal Enthalpies= -76.434769 Sum of electronic and thermal Free Energies= -76.456191 E (Thermal) CV S KCal/Mol Cal/Mol-Kelvin Cal/Mol-Kelvin Total 15.139 6.007 45.085 Electronic 0.000 0.000 0.000 Translational 0.889 2.981 34.608 Rotational 0.889 2.981 10.470 Vibrational 13.362 0.046 0.007 Q Log10(Q) Ln(Q) Total Bot 0.209715D-01 -1.678370 -3.864589 Total V=0 0.130223D+09 8.114687 18.684757 Vib (Bot) 0.161103D-09 -9.792897 -22.548979 Vib (V=0) 0.100037D+01 0.000160 0.000368 Electronic 0.100000D+01 0.000000 0.000000 Translational 0.300432D+07 6.477746 14.915562 Rotational 0.433292D+02 1.636781 3.768828 In Chem350 we went through all of the machinery that is involved to do the calculations. From the chem350 website you can download a Matlab file. In it you find the above data for the water molecule, and the calculation of thermodynamic properties. Going through the Matlab file and the Gaussian output file, you will find exact agreement between what we do in class, and what is implemented in the Gaussian program. These things are used in actual state-of-the-art calculations! I encourage you to go through the Matlab file and make sure you understand all the steps in the calculation. Find the corresponding results in the Gaussian output file.

6

ii) CO2 (linear molecule). To make sure we get things correctly also for linear molecules, let us do a calculation on the CO2 molecule, using the same DFT/B3LYP method using the cc-pVTZ basis set. This is the output I obtained, and need to set up the corresponding matlab calculation. From the Summary E(RB3LYP) -188.66056820 From the output file (view file): Frequencies -- 671.7935 671.7935 1371.8079 Frequencies -- 2417.0441 Molecular mass:

43.98983 amu.

Principal axes and moments of inertia in atomic units: 1 2 3 Eigenvalues -- 0.00000 153.80801 153.80801 Rotational symmetry number 2. This is all the information needed to calculate thermodynamical properties. You can feed it into the Matlab file and calculate properties directly. Then compare with the Gaussian thermochemistry output section: Zero-point correction= 0.011693 (Hartree/Particle) Thermal correction to Energy= 0.014311 Thermal correction to Enthalpy= 0.015255 Thermal correction to Gibbs Free Energy= -0.008999 Sum of electronic and zero-point Energies= -188.648876 Sum of electronic and thermal Energies= -188.646258 Sum of electronic and thermal Enthalpies= -188.645314 Sum of electronic and thermal Free Energies= -188.669567 E (Thermal) CV S KCal/Mol Cal/Mol-Kelvin Cal/Mol-Kelvin Total 8.980 6.855 51.046 Electronic 0.000 0.000 0.000 Translational 0.889 2.981 37.270 Rotational 0.592 1.987 13.073 Vibrational 7.499 1.887 0.703 Q Log10(Q) Ln(Q) Total Bot 0.137819D+05 4.139308 9.531109 Total V=0 0.329230D+10 9.517499 21.914852 Vib (Bot) 0.453971D-05 -5.342972 -12.302648 Vib (V=0) 0.108447D+01 0.035219 0.081095 Electronic 0.100000D+01 0.000000 0.000000 Translational 0.114679D+08 7.059484 16.255062 Rotational 0.264726D+03 2.422797 5.578695

7

D. Calculating equilibrium constants from Gibbs free energy calculations. Consider the reaction C2 H 2 ( g ) + H 2 ( g ) → C2 H 4 ( g ) To calculate the equilibrium constant at standard temperature and pressure one needs to 0 calculate the ! r G298.15 , i.e. the Gibbs free energy for each of the species. Such calculations are straightforward with what we know now. Create each of the molecular species in Gaussian and optimize geometries and calculate frequencies at for example the DFT/B3LYP/cc-pVTZ level. The most critical aspect of the calculation is the electronic energy. One can try to improve the accuracy in this respect. Since this is not a class on computational chemistry let us simply proceed at the above level. From the Gaussian outputs you can obtain the Gibbs free energy for each species, and also Cv and hence Cp using ideal gas expressions. From the calculated data calculate the equilibrium constant at 298.15 K, 1 atm, and also at 400 K, 1 atm. Also calculate the equilibrium constant using data in Reid and Engel (or from the web!). Compare the results. E. Calculation of transition states. Comparing Gaussian output with Matlab calculation. As a very simple, standard, example let us consider the hydgron migration reaction in HCN. You need to provide a guess for the transition state (elongate both CH and CN bonds a bit, in addition to changing the angle). Then select OPT+FREQ under jobtype, and select optimize to a transition state (Berny TS), and calculate the force constant matrix once. Under method select the DFT / B3LYP /cc-pVTZ methodology as usual. Upon running the calculation we can get the information we need. The first thing to do is check the Summay/Vibrations. We should find exactly one imaginary frequency (indicated in Gaussian as a negative frequency). You can look at it, and verify that this normal mode corresponds to the reaction coordinate. To do the statmech we lift the usual information from summary and output files: E(RB3LYP) -93.38546704 a.u. Frequencies -- -1129.0406 2066.9596 2592.7695 Molecular mass: 27.01090 amu. Principal axes and moments of inertia in atomic units: 1 2 3 Eigenvalues -- 4.38207 32.45085 36.83292 We can put this information in the Matlab file and run the thermochemistry. Note that the negative frequency does not enter the equations. Instead we get a factor kT/h in the rate constant. You can verify that the Gaussian output indeed agrees with our statmech calculation in Matlab. To calculate reaction rates and activation barriers we should also calculate the free energy of reactants (HCN) and products (HNC). This then provides the free energy reaction 8

profile, and a first estimate of the rate constant for the elementary reaction. Using these values the reaction rates can be calculated at any temperature, using the formulas for Gibbs free energy. You can compare with the traditional thermodynamics way of calculating the temperature and pressure dependence of the Gibbs free energy and the Arrhenius law. F. A case study: Comparing the rate of reaction for cis-trans isomerization in butadiene to the rate of conformational change in ethane: The effects of conjugation. The following calculations are time consuming!! Get them to run but you may wish to collect results later on. Since you can only run one Gaussian job it is all a bit painful. Such is life! Throughout we will use the DFT/B3LYP/cc-PVTZ approach. Optimize the structures of cis and trans butadiene and calculate vibrational frequencies. It is helpful to impose the correct pointgroup symmetry. This can be done using the pointgroup option under edit. Using the correct symmetry speeds up the calculation. If you want to change the structure and perhaps the symmetry, you have to undo setting the symmetry first. Also locate the transition state structure, by using an initial guess for the dihedral angel of 90 degrees and obtain vibrational frequencies. Check that you get one negative (=imaginary) frequency. Now you can process the free energies and calculate the forward and backward reaction rates, as well as the equilibrium constant at ambient temperature and pressure. You can also investigate the various contributions to the thermal properties (rotational, vibrational etc.). I anticipate the energy barrier to be substantial because the system is conjugated with its alternating pattern of double bonds. Let us compare to the barrier to rotation in ethane therefore. You can make the so-called staggered and eclipsed configurations of ethane, the latter being a transition state. Follow the same recipe and comment on the differences. See if you can find values for the above quantities in the literature or using the Web. For example visit the NIST side, which has a wealth of computational and experimental data. G. Summary All gas phase reactions can in principle be studied (given enough computer resources). We have now discussed all of the underlying theoretical principles in chem254 (thermodynamics), chem350 (kinetics and statmech) and chem356 (Quantum Mechanics): particle in the box: translational energy; Rigid Rotor model, (using spherical harmonics): rotational energy contributions; harmonic oscillator models for polyatomics: vibrational part of the energy. The major component of the reaction enthalpy is the difference in the electronic energy at the minimum of the potential (the optimized geometry), for each of the molecules involved in the reaction. Here accurate calculations are critical, and in prior classes we only scratched the surface regarding electronic energy calculations. You can obtain more background information on this in classes on computational chemistry (chem440 / computational), and quantum mechanics in

9

chemistry (also chem440 special topics), but it is a rather specialized subject without easy solutions.

10