Using the Output File from a Gaussian Frequency Calculation to Compute Ideal-gas Thermodynamic Functions

Ideal-gas thermodynamic functions such as molar entropy, heat capacity, and enthalpy content, can be computed easily from the molecular partition function. The details are described in many places, including here. The rigid-rotor/harmonic-oscillator (RRHO) approximation, adopted here, is popular because it is simple and is often adequate. All the necessary information is usually included in the output file from a vibrational frequency calculation executed using the quantum chemistry software from Gaussian, Inc. You can download a free Perl script (link at the bottom of this page), called thermo.pl, that will automatically extract the essential data from a Gaussian output file and compute thermodynamic functions at several temperatures. The basic data are also saved to a keyword-driven text file, which can be edited and used as input to the Perl script. 

This information can be used to compute quantities such as temperature-dependent free energies of reaction, delta-G(T).  Here is a worked example, for computing the relative stabilities of the cis- and trans- isomers of 1,2-dichloroethylene, CHCl=CHCl.  

---

If you don't want to install the script locally, it's now available on-line (as of Jan. 16, 2004).

---

Basic use of thermo.pl

Try the command "perl -v" on your computer. If it doesn't do anything, you will have to install Perl on your system before continuing.

If your frequency output file is called "chcl3_freq.out", then just type the following command:

perl thermo.pl chcl3_freq.out

You should get output like the following:

opening chcl3_freq.out for input

Scaling vibrational frequencies by 1.0000
Molecular mass = 117.914 u
External symmetry number = 3.
Rotational constants (GHz) = 3.19889 3.19889 1.65801 Harmonic vibrational frequencies (cm^-1):

	Unscaled	Scaled
1	261.0	261.0
2	261.0	261.0
3	367.4	367.4
4	668.1	668.1
5	737.9	737.9
6	737.9	737.9
7	1266.1	1266.1
8	1266.1	1266.1
9	3199.1	3199.1

Electronic States:
	degen	cm-1
1	1	0.0

Vibrational zero-point energy = 4382.2 cm^-1 = 52.42 kJ/mol.

T(K)	S (J/mol.K)	Cp (J/mol.K)	ddH (kJ/mol)
100.00	239.44	40.30	3.50
200.00	271.96	54.64	8.27
298.15	295.94	65.76	14.20
300.00	296.35	65.94	14.33
400.00	316.55	74.43	21.37
500.00	333.85	80.61	29.13
600.00	348.97	85.14	37.43
700.00	362.36	88.56	46.13
800.00	374.37	91.24	55.12
900.00	385.25	93.40	64.36
1000.00	395.18	95.18	73.79

Input data (including scaled vibrational frequencies) written to file 'thermo.dat'.

The column headings are:

T (temperature in K)
S (entropy in J mol-1 K-1)
Cp (heat capacity at constant pressure in kJ mol-1 K-1)
ddH = H(T) - H(0) (enthalpy content, in kJ mol-1).

For comparison, here are the corresponding values from the JANAF Tables (by M. W. Chase; 4th edition): 

T(K)	S (J/mol.K)	Cp (J/mol.K)	ddH (kJ/mol)
100.00	239.35	40.32	3.50
200.00	271.80	54.30	8.26
298.15	295.62	65.38	14.15
300.00	296.03	65.57	14.27
400.00	316.14	74.26	21.29
500.00	333.43	80.61	29.05
600.00	348.56	85.26	37.35
700.00	361.98	88.77	46.06
800.00	374.02	91.51	55.08
900.00	384.93	93.72	64.35
1000.00	394.90	95.52	73.81

If the file "thermo.dat" already exists, you will be prompted to overwrite it or not; the default is not to overwrite. "Thermo.dat" contains the important information and can be used as input to the script:

perl thermo.pl thermo.dat

This command should yield the same results as before.

If you want a temperature (e.g., 550 K) that is not in the default list, supply that temperature as an additional command-line argument:

perl thermo.pl chcl3_freq.out 550

Likewise, if you want to scale the vibrational frequencies by a multiplicative factor (e.g., 1.0015), specify that factor as the additional command-line argument, as shown below. The program examines the magnitude of the numerical arguments to distinguish a temperature (greater than 10) from a vibrational scaling factor (between 0.5 and 1.5).

perl thermo.pl chcl3_freq.out 1.0015

You can specify both a scaling factor and an extra temperature:

perl thermo.pl chcl3_freq.out 550 1.0015

Advanced use: editing data files

Additional functions of thermo.pl are only available using a keyword-driven data file as input. The data file must be edited to supply additional keywords and the corresponding data.

In the example above, the program automatically generates a data file (thermo.dat) that looks like this:

MASS	117.914380
GHZ	3.198890	3.198890	1.658010
VIB	261.0	261.0	367.4	668.1	737.9	737.9	1266.1	1266.1	1266.1
SYMNO	3
ELEC	1	0.0

Here are the keywords that the program understands. Each may be repeated as many times as makes sense. Data files can be created in any way; no Gaussian output files are actually necessary for using thermo.pl. 

Keyword	Meaning	Unit	Data required
MASS	molecular mass	atomic units (u)	mass
GHZ	rotational constant(s)	GHz (1 cm-1 ~ 29.979 GHz)	one constant for a linear molecule, three constants for a nonlinear molecule
VIB	harmonic vibrational frequency	cm-1	list of frequencies
SYNMO	symmetry number	none	symmetry number
ELEC	electronic energy level	none for 1st argument; cm-1 for 2nd	degeneracy and energy
AVID	anharmonic vibrational ladder	cm-1	harmonic constant, anharmonicity constant, and optional high-energy limit (for truncation)
XVIB	explicit energy levels	cm-1	list of energy values; the lowest level should be at 0 energy

Example: HF molecule

MASS	20.006230
GHZ	628.23
AVIB	4138.73	90.05
SYMNO	1
ELEC	1	0.0

The program will apply a reasonable limit on anharmonic vibrational ladders, even if none is specified (as here). At 298.15 K, here are the results from the above file and from the 4th edition of the JANAF tables:

Source	S (J/mol.K)	Cp (J/mol.K)	ddH (kJ/mol)
thermo.pl	173.59	29.10	8.68
JANAF4	173.78	29.14	8.60

Example: F atom

A Gaussian frequency calculation on the F atom will produce an output file suitable for thermo.pl. However, the script only detects spin degeneracy, not spatial degeneracy. So the degeneracy of the (ground) electronic level will come out as 2, which is wrong. If spin-orbit coupling is ignored, the degeneracy should be 6 for the ²P ground state. Of course, the spin-orbit splitting is important (404.1 cm^-1) and should not be ignored. Here are the data file and the comparison with JANAF:

MASS	18.998400
GHZ
SYMNO	1
ELEC	4	0.0
ELEC	2	404.1

Source	S (J/mol.K)	Cp (J/mol.K)	ddH (kJ/mol)
thermo.pl	158.75	22.75	6.52
JANAF4	158.75	22.75	6.52

Downloading thermo.pl

You may download and use the program without charge or registration. Do not distribute the software without attribution and this documentation.

To download. 

Citation

K.K. Irikura, THERMO.PL, National Institute of Standards and Technology, 2002.