Calculating free energies of constrained species in Gaussian

scripts to perform frequency analysis and compute free energy of surface species (species with constrained atoms) in Gaussian

Introduction

Gaussian09 has inbuilt features to calculate the free energies (translational, rotational, and vibrational) of gas phase species. However, it does not readily calculate the free energies of species with constrained atoms. An example of this can be when we are trying to calculate the free energy of surface species. The surface is at times represented by a small cluster, where the cluster is capped by terminating atoms and some peripheral atoms of the cluster are held fixed to mimic the rigidity of the surface. The figure below shows a vicinal silanol cluster used to represent a vicinal silanol site on the surface of silica.

Here, the cluster is capped by hydrogen atoms. And the peripheral OH groups (shown in red) are held fixed. To calculate the vibrational free energy we expand the enregy around the minimum (structure corresponding to a local minima) as follows:

Here, ∆x is the deviation from the minima, H is the hessian computed at the minima, and E₀ is the energy of the minima. H is written as follows:

.

Here, x_i is the ith coordinate. A molecule with N atoms has a total of 3N coordinates. Now, if we divide the hessian matrix into coordinates of peripheral (fixed) and non-peripheral (free) atoms, we can represent it as follows:

.

Here, the red portion refers to peripheral atoms (fixed) and blue section refers to non-peripheral atoms (free). To eliminate the effects of the peripheral atoms only the blue sub-matrix is considered and is hence referred to as the reduced hessian represented as H_red.

The reduced hessian is mass weighted and diagonalized as follows:

.

Here, H_MW is the reduced mass-weighted hessian and M is the mass matrix represented as

.

Here, M_i is the mass of the atom corresponding to coordinate x_i. Following this the normal vibrational frequencies can be calculated using the eigenvalues of H_MW

.

Here ν_i is the freq1uency of vibration of the ith normal mode and ε_i is the ith eigenvalue of H_MW

The vibrational frequencies can be used to compute the ith vibrational partition function as follows:

.

The total vibrational partition function is written as the product of the individual partition functions

.

The vibrational partition function can then be used to calculate the vibrational contribution to the free energy

.

Usage

\g_free_e contains the code and an example Gaussian output file (vicinal silica site in Fig. 1).

Place all the atoms to be fixed at the end in the Guassian input file (opt.gjf).

When running a gaussian calculation generate a checkpoint file and convert it to a formatted checkpoint file using the following:

<checkpoint filename> formchk <formatted checkpoint filename>

The formatted checkpoint file contains the second order derivatives (Hessian). Following this run the generete_report.py as follows:

python .\free_e_module.py <number of fixed atoms> <temperature> <output gaussian filename> <report filename>

For example,

python .\free_e_module.py 8 298.15 opt.log Report.txt

The output file reports to total free energy, enthalpy, entropy and normal vibrational frequencies.

References

[1] F. Jensen. Introduction to Computational Chemistry. Wiley (1999)