## Thermochemistry in Gaussian

Joseph W. Ochterski, Ph.D.
help@gaussian.com

April 19, 2000

Abstract:

The purpose of this paper is to explain how various thermochemical values are computed in Gaussian. The paper documents what equations are used to calculate the quantities, but doesn't explain them in great detail, so a basic understanding of statistical mechanics concepts, such as partition functions, is assumed. Gaussian thermochemistry output is explained, and a couple of examples, including calculating the enthalpy and Gibbs free energy for a reaction, the heat of formation of a molecule and absolute rates of reaction are worked out.

Contents

Introduction

The equations used for computing thermochemical data in Gaussian are equivalent to those given in standard texts on thermodynamics. Much of what is discussed below is covered in detail in ``Molecular Thermodynamics'' by McQuarrie and Simon (1999). I've cross-referenced several of the equations in this paper with the same equations in the book, to make it easier to determine what assumptions were made in deriving each equation. These cross-references have the form [McQuarrie, §7-6, Eq. 7.27] which refers to equation 7.27 in section 7-6.

One of the most important approximations to be aware of throughout this analysis is that all the equations assume non-interacting particles and therefore apply only to an ideal gas. This limitation will introduce some error, depending on the extent that any system being studied is non-ideal. Further, for the electronic contributions, it is assumed that the first and higher excited states are entirely inaccessible. This approximation is generally not troublesome, but can introduce some error for systems with low lying electronic excited states.

The examples in this paper are typically carried out at the HF/STO-3G level of theory. The intent is to provide illustrative examples, rather than research grade results.

The first section of the paper is this introduction. The next section of the paper, I give the equations used to calculate the contributions from translational motion, electronic motion, rotational motion and vibrational motion. Then I describe a sample output in the third section, to show how each section relates to the equations. The fourth section consists of several worked out examples, where I calculate the heat of reaction and Gibbs free energy of reaction for a simple bimolecular reaction, and absoloute reaction rates for another. Finally, an appendix gives a list of the all symbols used, their meanings and values for constants I've used.

Sources of components for thermodynamic quantities

In each of the next four subsections of this paper, I will give the equations used to calculate the contributions to entropy, energy, and heat capacity resulting from translational, electronic, rotational and vibrational motion. The starting point in each case is the partition function q(V,T) for the corresponding component of the total partition function. In this section, I'll give an overview of how entropy, energy, and heat capacity are calculated from the partition function.

The partition function from any component can be used to determine the entropy contribution S from that component, using the relation [McQuarrie, §7-6, Eq. 7.27]:

The form used in Gaussian is a special case. First, molar values are given, so we can divide by , and substitute . We can also move the first term into the logarithm (as e), which leaves (with N=1):

The internal thermal energy E can also be obtained from the partition function [McQuarrie, §3-8, Eq. 3.41]:

and ultimately, the energy can be used to obtain the heat capacity [McQuarrie, §3.4, Eq. 3.25]:

These three equations will be used to derive the final expressions used to calculate the different components of the thermodynamic quantities printed out by Gaussian.

Contributions from translation

The equation given in McQuarrie and other texts for the translational partition function is [McQuarrie, §4-1, Eq. 4.6]:

The partial derivative of with respect to T is:

which will be used to calculate both the internal energy and the third term in Equation 1.

The second term in Equation 1 is a little trickier, since we don't know V. However, for an ideal gas, , and . Therefore,

which is what is used to calculate in Gaussian. Note that we didn't have to make this substitution to derive the third term, since the partial derivative has V held constant.

The translational partition function is used to calculate the translational entropy (which includes the factor of e which comes from Stirling's approximation):

The contribution to the internal thermal energy due to translation is:

Finally, the constant volume heat capacity is given by:

Contributions from electronic motion

The usual electronic partition function is [McQuarrie, §4-2, Eq. 4.8]:

where is the degeneracy of the energy level, is the energy of the n-th level.

Gaussian assumes that the first electronic excitation energy is much greater than . Therefore, the first and higher excited states are assumed to be inaccessible at any temperature. Further, the energy of the ground state is set to zero. These assumptions simplify the electronic partition function to:

which is simply the electronic spin multiplicity of the molecule.

The entropy due to electronic motion is:

Since there are no temperature dependent terms in the partition function, the electronic heat capacity and the internal thermal energy due to electronic motion are both zero.

Contributions from rotational motion

The discussion for molecular rotation can be divided into several cases: single atoms, linear polyatomic molecules, and general non-linear polyatomic molecules. I'll cover each in order.

For a single atom, . Since does not depend on temperature, the contribution of rotation to the internal thermal energy, its contribution to the heat capacity and its contribution to the entropy are all identically zero.

For a linear molecule, the rotational partition function is [McQuarrie, §4-6, Eq. 4.38]:

where . I is the moment of inertia. The rotational contribution to the entropy is

The contribution of rotation to the internal thermal energy is

and the contribution to the heat capacity is

For the general case for a nonlinear polyatomic molecule, the rotational partition function is [McQuarrie, §4-8, Eq. 4.56]:

Now we have , so the entropy for this partition function is

Finally, the contribution to the internal thermal energy is

and the contribution to the heat capacity is

The average contribution to the internal thermal energy from each rotational degree of freedom is RT/2, while it's contribution to is R/2.

Contributions from vibrational motion

The contributions to the partition function, entropy, internal energy and constant volume heat capacity from vibrational motions are composed of a sum (or product) of the contributions from each vibrational mode, K. Only the real modes are considered; modes with imaginary frequencies (i.e. those flagged with a minus sign in the output) are ignored. Each of the (or for linear molecules) modes has a characteristic vibrational temperature, .

There are two ways to calculate the partition function, depending on where you choose the zero of energy to be: either the bottom of the internuclear potential energy well, or the first vibrational level. Which you choose depends on whether or not the contributions arising from zero-point energy will be coputed separately or not. If they are computed separately, then you should use the bottom of the well as your reference point, otherwise the first vibrational energy level is the appropriate choice.

If you choose the zero reference point to be the bottom of the well ( BOT), then the contribution to the partition function from a given vibrational mode is [McQuarrie, §4-4, Eq. 4.24]:

and the overall vibrational partition function is[McQuarrie, §4-7, Eq. 4.46]:

On the other hand, if you choose the first vibrational energy level to be the zero of energy (V=0), then the partition function for each vibrational level is

and the overall vibrational partition function is:

Gaussian uses the bottom of the well as the zero of energy (BOT) to determine the other thermodynamic quantities, but also prints out the V=0 partition function. Ultimately, the only difference between the two references is the additional factor of , (which is the zero point vibrational energy) in the equation for the internal energy . In the expressions for heat capacity and entropy, this factor disappears since you differentiate with respect to temperature (T).

The total entropy contribution from the vibrational partition function is:

To get from the fourth line to the fifth line in the equation above, you have to multiply by .

The contribution to internal thermal energy resulting from molecular vibration is

Finally, the contribution to constant volume heat capacity is

Low frequency modes (defined below) are included in the computations described above. Some of these modes may be internal rotations, and so may need to be treated separately, depending on the temperatures and barriers involved. In order to make it easier to correct for these modes, their contributions are printed out separately, so that they may be subtracted out. A low frequency mode in Gaussian is defined as one for which more than five percent of an assembly of molecules are likely to exist in excited vibrational states at room temperature. In other units, this corresponds to about 625 cm ,  Hz, or a vibrational temperature of 900 K.

It is possible to use Gaussian to automatically perform some of this analysis for you, via the Freq=HindRot keyword. This part of the code is still undergoing some improvements, so I will not go into it in detail. See P. Y. Ayala and H. B. Schlegel, J. Chem. Phys. 108 2314 (1998) and references therein for more detail about the hindered rotor analysis in Gaussian and methods for correcting the partition functions due to these effects.

## Thermochemistry output from Gaussian

This section describes most of the Gaussian thermochemistry output, and how it relates to the equations I've given above.

Output from a frequency calculation

In this section, I intentionally used a non-optimized structure, to show more output. For production runs, it is very important to use structures for which the first derivatives are zero -- in other words, for minima, transition states and higher order saddle points. It is occasionally possible to use structures where one of the modes has non-zero first derivatives, such as along an IRC. For more information about why it is important to be at a stationary point on the potential energy surface, see my white paper on ``Vibrational Analysis in Gaussian ''.

Much of the output is self-explanatory. I'll only comment on some of the output which may not be immediately clear. Some of the output is also described in ``Exploring Chemistry with Electronic Structure Methods, Second Edition'' by James B. Foresman and Æleen Frisch.

``` -------------------
- Thermochemistry -
-------------------
Temperature   298.150 Kelvin.  Pressure   1.00000 Atm.
Atom  1 has atomic number  6 and mass  12.00000
Atom  2 has atomic number  6 and mass  12.00000
Atom  3 has atomic number  1 and mass   1.00783
Atom  4 has atomic number  1 and mass   1.00783
Atom  5 has atomic number  1 and mass   1.00783
Atom  6 has atomic number  1 and mass   1.00783
Atom  7 has atomic number  1 and mass   1.00783
Atom  8 has atomic number  1 and mass   1.00783
Molecular mass:    30.04695 amu.```

The next section gives some characteristics of the molecule based on the moments of inertia, including the rotational temperature and constants. The zero-point energy is calculated using only the non-imaginary frequencies.

``` Principal axes and moments of inertia in atomic units:
1           2           3
EIGENVALUES --    23.57594    88.34097    88.34208
X            1.00000     0.00000     0.00000
Y            0.00000     1.00000     0.00001
Z            0.00000    -0.00001     1.00000
THIS MOLECULE IS AN ASYMMETRIC TOP.
ROTATIONAL SYMMETRY NUMBER    1.
ROTATIONAL TEMPERATURES (KELVIN)     3.67381     0.98044    0.98043
ROTATIONAL CONSTANTS (GHZ)          76.55013    20.42926   20.42901
Zero-point vibrational energy    204885.0 (Joules/Mol)
48.96870 (Kcal/Mol)```

If you see the following warning, it can be a sign that one of two things is happening. First, it often shows up if your structure is not a minimum with respect to all non-imaginary modes. You should go back and re-optimize your structure, since all the thermochemistry based on this structure is likely to be wrong. Second, it may indicate that there are internal rotations in your system. You should correct for errors caused by this situation.

``` WARNING-- EXPLICIT CONSIDERATION OF    1 DEGREES OF FREEDOM AS
VIBRATIONS MAY CAUSE SIGNIFICANT ERROR```

Then the vibrational temperatures and zero-point energy (ZPE):

``` VIBRATIONAL TEMPERATURES:     602.31   1607.07   1607.45    1683.83   1978.85
(KELVIN)          1978.87   2303.03   2389.95   2389.96   2404.55
2417.29   2417.30   4202.52   4227.44   4244.32
4244.93   4291.74   4292.31

Zero-point correction=                      0.078037 (Hartree/Particle)```

Each of the next few lines warrants some explanation. All of them include the zero-point energy. The first line gives the correction to the internal thermal energy, .

` Thermal correction to Energy=             0.081258`

The next two lines, respectively, are

and

where .

``` Thermal correction to Enthalpy=                 0.082202
Thermal correction to Gibbs Free Energy=        0.055064```

The Gibbs free energy includes , so when it's applied to calculating for a reaction, is already included. This means that will be computed correctly when the number of moles of gas changes during the course of a reaction.

The next four lines are estimates of the total energy of the molecule, after various corrections are applied. Since I've already used E to represent internal thermal energy, I'll use for the total electronic energy.

Sum of electronic and zero-point energies=

Sum of electronic and thermal energies=

Sum of electronic and thermal enthalpies=

Sum of electronic and thermal free energies=

``` Sum of electronic and zero-point Energies=           -79.140431
Sum of electronic and thermal Energies=              -79.137210
Sum of electronic and thermal Enthalpies=            -79.136266
Sum of electronic and thermal Free Energies=         -79.163404```

The next section is a table listing the individual contributions to the internal thermal energy ( ), constant volume heat capacity ( ) and entropy ( ). For every low frequency mode, there will be a line similar to the last one in this table (labeled VIBRATION 1). That line gives the contribution of that particular mode to , and . This allows you to subtract out these values if you believe they are a source of error.

```                     E (Thermal)             CV                S
KCAL/MOL        CAL/MOL-KELVIN    CAL/MOL-KELVIN
TOTAL                   50.990              8.636             57.118
ELECTRONIC               0.000              0.000              0.000
TRANSLATIONAL            0.889              2.981             36.134
ROTATIONAL               0.889              2.981             19.848
VIBRATIONAL             49.213              2.674              1.136
VIBRATION  1             0.781              1.430              0.897```

Finally, there is a table listing the individual contributions to the partition function. The lines labeled BOT are for the vibrational partition function computed with the zero of energy being the bottom of the well, while those labeled with (V=0) are computed with the zero of energy being the first vibrational level. Again, special lines are printed out for the low frequency modes.

```                       Q             LOG10(Q)            LN(Q)
TOTAL BOT       0.470577D-25       -25.327369       -58.318422
TOTAL V=0       0.368746D+11        10.566728        24.330790
VIB (BOT)       0.149700D-35       -35.824779       -82.489602
VIB (BOT)  1    0.419879D+00        -0.376876        -0.867789
VIB (V=0)       0.117305D+01         0.069318         0.159610
VIB (V=0)  1    0.115292D+01         0.061797         0.142294
ELECTRONIC      0.100000D+01         0.000000         0.000000
TRANSLATIONAL   0.647383D+07         6.811161        15.683278
ROTATIONAL      0.485567D+04         3.686249         8.487901```

Output from compound model chemistries

This section explains what the various thermochemical quantities in the summary out of a compound model chemistry, such as CBS-QB3 or G2, means.

I'll use a CBS-QB3 calculation on water, but the discussion is directly applicable to all the other compound models available in Gaussian. The two lines of interest from the output look like:

``` CBS-QB3 (0 K)=           -76.337451 CBS-QB3 Energy=             -76.334615
CBS-QB3 Enthalpy=        -76.333671 CBS-QB3 Free Energy=        -76.355097```

Here are the meanings of each of those quantities.

• CBS-QB3 (0 K): This is the total electronic energy as defined by the compound model, including the zero-point energy scaled by the factor defined in the model.
• CBS-QB3 Energy: This is the total electronic energy plus the internal thermal energy, which comes directly from the thermochemistry output in the frequency part of the output. The scaled zero-point energy is included.
• CBS-QB3 Enthalpy: This is the total electronic energy plus , as it is described above with the unscaled zero-point energy removed. This sum is appropriate for calculating enthalpies of reaction, as described below.
• CBS-QB3 Free Energy: This is the total electronic energy plus , as it is described above, again with the unscaled zero-point energy removed. This sum is appropriate for calculating Gibbs free energies of reaction, also as described below.

## Worked-out Examples

In this section I will show how to use these results to generate various thermochemical information.

I've run calculations for each of the reactants and products in the reaction where ethyl radical abstracts a hydrogen atom from molecular hydrogen:

as well as for the transition state (all at 1.0 atmospheres and 298.15K). The thermochemistry output from Gaussian is summarized in Table 1.

Table 1: Calculated thermochemistry values from Gaussian for the reaction . All values are in Hartrees.

Once you have the data for all the relevant species, you can calculate the quantities you are interested in. Unless otherwise specified, all enthalpies are at 298.15K. I'll use to shorten the equations.

Enthalpies and Free Energies of Reaction

The usual way to calculate enthalpies of reaction is to calculate heats of formation, and take the appropriate sums and difference.

However, since Gaussian provides the sum of electronic and thermal enthalpies, there is a short cut: namely, to simply take the difference of the sums of these values for the reactants and the products. This works since the number of atoms of each element is the same on both sides of the reaction, therefore all the atomic information cancels out, and you need only the molecular data.

For example, using the information in Table 1, the enthalpy of reaction can be calculated simply by

The same short cut can be used to calculate Gibbs free energies of reaction:

Rates of Reaction

In this section I'll show how to compute rates of reaction using the output from Gaussian. I'll be using results derived from transition state theory in section 28-8 of ``Physical Chemistry, A Molecular Approach'' by D. A. McQuarrie and J. D. Simon. The key equation (number 28.72, in that text) for calculating reaction rates is

I'll use for the concentration. For simple reactions, the rest is simply plugging in the numbers. First, of course, we need to get the numbers. I've run HF/STO-3G frequency calculation for reaction FH + Cl F + HCl and also for the reaction with deuterium substituted for hydrogen. The results are summarized in Table 2. I put the total electronic energies of each compound into the table as well to illustrate the point that the final geometry and electronic energy are independent of the masses of the atoms. Indeed, the cartesian force constants themselves are independent of the masses. Only the vibrational analysis, and quantities derived from it, are mass dependent.

Table 2: Total HF/STO-3G electronic energies and sum of electronic and thermal free energies for atoms molecules and transition state complex of the reaction FH + Cl F + HCl and the deuterium substituted analog.

The first step in calculating the rates of these reactions is to compute the free energy of activation, ((H) is for the hydrogen reaction, (D) is for deuterium reaction):

Then we can calculate the reaction rates. The values for the constants are listed in the appendix. I've taken .

So we see that the deuterium reaction is indeed slower, as we would expect. Again, these calculations were carried out at the HF/STO-3G level, for illustration purposes, not for research grade results. More complex reactions will need more sophisticated analyses, perhaps including careful determination of the effects of low frequency modes on the transition state, and tunneling effects.

Enthalpies and Free Energies of Formation

Calculating enthalpies of formation is a straight-forward, albeit somewhat tedious task, which can be split into a couple of steps. The first step is to calculate the enthalpies of formation ( ) of the species involved in the reaction. The second step is to calculate the enthalpies of formation of the species at 298K.

Calculating the Gibbs free energy of reaction is similar, except we have to add in the entropy term:

To calculate these quantities, we need a few component pieces first. In the descriptions below, I will use M to stand for the molecule, and X to represent each element which makes up M, and x will be the number of atoms of X in M.

• Atomization energy of the molecule, :

These are readily calculated from the total energies of the molecule ( ), the zero-point energy of the molecule ( ) and the constituent atoms:

• Heats of formation of the atoms at 0K, :

I have tabulated recommended values for the heats of formation of the first and second row atomic elements at 0K in Table 3. There are (at least) two schools of thought with respect to the which atomic heat of formation data is most appropriate. Some authors prefer to use purely experimental data for the heats of formation (Curtiss, et. al., J.\ Chem. Phys. 106, 1063 (1997)). Others (Ochterski, et. al., J.\ Am. Chem. Soc. 117, 11299 (1995)), prefer to use a combination of experiment and computation to obtain more accurate values. The elements of concern are boron, beryllium and silicon. I have arranged Table 3 so that either may be used; the experimental results are on the top, while the computational results for the three elements with large experimental uncertainty are listed below that.

• Enthalpy corrections of the atomic elements, :

The enthalpy corrections of the first and second row atomic elements are also included in Table 3. These are used to convert the atomic heats of formation at 0K to those at 298.15K, and are given for the elements in their standard states. Since they do not depend on the accuracy of the heat of formation of the atom, they are the same for both the calculated and experimental data.

This is generally not the same as the output from Gaussian for a calculation on an isolated gas-phase atom. These values are referenced to the standard states of the elements. For example, the value for hydrogen atom is 1.01 kcal/mol. This is , and not , which is what Gaussian calculates.

Table: Experimental and calculated enthalpies of formation of elements (kcal mol ) and entropies of atoms (cal mol K ) .  Experimental enthalpy values taken from J. Chem. Phys. 106, 1063 (1997).  Calculated enthalpy values taken from J. Am.\ Chem. Soc. 117, 11299 (1995). is the same for both calculated and experimental, and is taken for elements in their standard states. Entropy values taken from JANAF Thermochemical Tables: M. W. Chase, Jr., C. A. Davies, J. R. Downey, Jr., D. J. Frurip, R. A. McDonald, and A. N. Syverud, J. Phys. Ref. Data 14 Suppl. No. 1 (1985). No error bars are given for Mg in JANAF.

• Enthalpy correction for the molecule, :

For a molecule, this is , where is the value printed out in the line labeled `` Thermal correction to Enthalpy'' in Gaussian output. Remember, though, that the output is in Hartrees/particle, and needs to be converted to kcal/mol. The conversion factor is 1 Hartree=627.5095 kcal/mol.

• Entropy for the atoms, :

These are summarized in Table 3. The values were taken from the JANAF tables (M. W. Chase, Jr., C. A. Davies, J. R. Downey, Jr., D. J. Frurip, R. A. McDonald, and A. N. Syverud, J. Phys. Ref. Data 14 Suppl. No. 1 (1985)). Again, these values are referenced to the standard states of the elements, and will disagree with values Gaussian calculates for the gas-phase isolated atoms.

• Entropy for the molecule, :

Gaussian calculates G=H - TS, and prints it out on the line labeled ``Thermal correction to Gibbs Free Energy''. The entropy can be teased out of this using S=(H - G)/T.

Putting all these pieces together, we can finally take the steps necessary to calculate and :

1. Calculate for each molecule:

2. Calculate for each molecule:

3. Calculate for each molecule:

Here is a worked out example, where I've calculated for the reactants and products in the reaction I described above.

First, I'll calculate for one of the species:

The next step is to calculate the :

To calculate the Gibbs free energy of formation, we have (the factors of 1000 are to convert kcal to or from cal) :

## Summary

The essential message is this: the basic equations used to calculate thermochemical quantities in Gaussian are based on those used in standard texts. Since the vibrational partition function depends on the frequencies, you must use a structure that is either a minimum or a saddle point. For electronic contributions to the partition function, it is assumed that the first and all higher states are inaccessible at the temperature the calculation is done at. The data generated by Gaussian can be used to calculate heats and free energies of reactions as well as absolute rate information.

Appendix: Symbols

=contribution to the heat capacity due to electronic motion
=contribution to the heat capacity due to rotational motion
=total heat capacity ( )
=contribution to the heat capacity due to translation
=contribution to the heat capacity due to vibrational motion
=internal energy due to electronic motion
=internal energy due to rotational motion
=total internal energy ( )
=internal energy due to translation
=internal energy due to vibrational motion
=correction to the Gibbs free energy due to internal energy
=correction to the enthalpy due to internal energy
I
=moment of inertia
K
=index of vibrational modes
N
=number of moles
P
=pressure (default is 1 atmosphere)
R
=gas constant
=entropy due to electronic motion
=entropy due to rotational motion
=total entropy ( )
=entropy due to translation
=entropy due to vibrational motion
T
=temperature (default is 298.15)
V
=volume
=characteristic temperature for rotation (in the x, y or z plane)
=characteristic temperature for vibration K
=Boltzmann constant
=the energy of the n-th energy level
=the degeneracy of the n-th energy level
=symmetry number for rotation
h
=Planck's constant
m
=mass of the molecule
n
=number of particles (always 1)
=electronic partition function
=rotational partition function
=translational partition function
=vibrational partition function
=zero-point energy of the molecule
=the total electronic energy, the MP2 energy for example
=standard enthalpy, entropy and Gibbs free energy -- each compound is in its standard state at the given temperature
=Gibbs free energy of activation
=concentration (taken to be 1)
k(T)
=reaction rate at temperature T