Selecting a Structure
One of the first conditions that must be satisfied in order to perform an accurate calculation of vibrational frequencies is that the underlying molecular structure represents a stationary point, that is, all first derivatives of the energy with respect to nuclear coordinates must be essentially zero. Using the vibrational frequency spectrum for the dihydrogen molecule (H-H) as an example, one can show that small deviations from the equilibrium structure can quickly lead to large errors in the calculated vibrational frequencies. The experimentally measured fundamental vibrational mode is located at 4160 cm^-1. Expanding the vibrational frequencies in a series containing harmonic and anharmonic terms, the corresponding (experimental) harmonic frequency is 4401 cm^-1. The following values are obtained at the HF/6-31G level of theory using structures optimized to different accuracy:

E (Hartree) r(H-H) (pm) rms gradient (Hartree/Bohr) frequencies (cm-1) -1.12680935484 73.500 3.86*10-3 -0.0001, -0.0001, 0.0001, 381.6219, 381.6219, 4570.76 -1.12682477931 73.200 1.16*10-3 -0.0002, -0.0001, -0.0001, 244.4793, 244.4793, 4615.51 -1.12682703185 73.100 8.06*10-4 -0.0001, 0.0001, 0.0002, 174.9713, 174.9713, 4630.52 -1.12682761115 73.050 4.20*10-4 -0.0001, -0.0001, 0.0000, 126.2380, 126.2380, 4638.04 -1.12682782422 73.000 3.13*10-5 -0.0001, 0.0000, 0.0001, 34.4845, 34.4845, 4645.57 -1.12682782541 72.996 1.86*10-7 -0.0002, 0.0000, 0.0001, 2.6607, 2.6607, 4646.2

All six structures summarized in the table above have rather similar total energies, the energy difference between the most favorable structure (last entry) and the least favorable structure (first entry) being less than 1 kJ/mol. The structural differences for these two structures amount to less than 1% (only 0.5 pm) in the H-H bond distance, giving rise to a residual gradient for the least favorable structure of 3.86*10^-3au. Compared to these differences in structure, energy and gradient, the changes in the calculated vibrational frequencies are much more pronounced. For a linear molecule with two atoms we expect only one real vibrational mode in addition to five near-zero modes describing translational and rotational motion. This latter expectation is not fullfilled for the first five entries in the table as in particular the two rotational modes are far from zero. The H-H stretching mode also suffers from the residual gradient, showing variations of almost 80 cm^-1. That the predicted harmonic frequency for the H-H streching mode appears to converge to a value further away from the experimentally measured one with decreasing residual gradient is primarily due to deficiencies of both the Hartree-Fock model as well as the small basis set (see below), but not a consequence of insufficient geometry optimization.

Similar results are obtained with other theoretical methods whose second derivatives cannot be calculated analytically and thus require a numerical algorithm. This is the case for the semiempirical AM1 method:

E (Hartree) r(H-H) (pm) rms gradient (Hartree/Bohr) frequencies (cm-1) -0.0081876 66.60 7.34*10-3 -550.9, -550.9, -0.0001, 0.0001, 0.0001, 4428.7 -0.0082603 67.60 4.01*10-4 -119.9, -119.9, -0.0001, 0.0001, 0.0001, 4346.4 -0.0082606 67.65 6.08*10-5 -20.8, -20.8, -0.0002, -0.0001, 0.0001, 4342.3 -0.0082606 67.66 1.01*10-6 -0.0001, 0.000, 0.0001, 44.907, 44.907, 4341.5

We observe that the AM1 value converges to 4342 cm^-1 with increasing accuracy of the geometry optimization. Again we note a large variation in predicted vibrational frequencies with an increasing residual gradient. In contrast to the HF results discussed before the five translational and rotational modes cannot be converged to very small values now. This is a consequence of additional problems of numerical vibrational frequency calculations (see below).

Selecting an Algorithm
Calculation of the force constant matrix (Hessian) can be performed with two different algorithms:

1) Using an analytical algorithm. The keyword for this option in Gaussian is

freq=Analytic

and it is the default for those theoretical methods, for which analytical second derivatives are available. The latter include Hartree-Fock theory (HF), all density functional methods (DFT), as well as MP2. Analytical calculation of the force constants is significantly faster and also more accurate than the alternative approach:

2) Using a numerical method. The keyword for this option in Gaussian is

freq=Numerical

and it is the default for methods for which analytical second derivatives are not available in Gaussian such as all semiempirical methods. For this option, the algorithm involves double displacements of all centers in all three cartesian coordinates. The step size used in this procedure can be specified with

freq=(Numerical,Step=N)

N being the size of the step in units of 0.0001 Angstroms. The default value for N is 100 in semiempirical calclations. Smaller values require a somewhat tighter convergence criterion for the energy calculations as the unperturbed and perturbed structures are structurally and energetically rather similar. This can be specified with

scf=(Conver=n)

which sets the convergence criterion in energy calculations to 10^-n Hartree on two successive SCF iterations. The default for semiempirical calculations such as AM1 is n=6, but values of n=7 or 8 are more appropriate once the step size in numerical second derivative calculations is less than N=50. Using these options, the vibrational frequency calculation for H₂ (AM1 level, H-H distance = 67.659 pm) can be modified as follows (the principal stretching vibration is again shown in bold):

keywords	frequencies (cm-1)
freq=(numerical,step=100) scf=(conver=6)	-0.0001, -0.0001, 0.0001, 45.2207, 45.2207, 4341.5252
freq=(numerical,step=50) scf=(conver=6)	-0.0001, 0.0001, 0.0001, 22.3738, 22.3738, 4341.2912
freq=(numerical,step=10) scf=(conver=7)	-0.0001, -0.0001, -0.0001, 1.6958, 1.6958, 4341.2165

The last entry in this table marks the limits of numerical accuracy that can be achieved with this approach.

Differentiating Ground and Transiton States
Pushing the limits of vibrational frequency calculations is sometimes necessary in order to achieve unequivocal characterization of minima and transition states on the potential energy surface. While minima are expected to show 6 very low frequencies in their vibrational frequency spectrum (3 for translational and 3 for rotational motions in a non-linear system) one large negative and 6 very small frequencies are expected for a transition state. The vibrational spectra calculated for methanol (CH₃OH) will be used as an example. For this system one would expect 3*6 - 6 = 12 true vibrational modes. The lowest seven vibrational modes are shown in the following tables, first for the staggered and then for the eclipsed conformer (HF/6-31G(d) results):

E (Hartree) rms gradient (Hartree/Bohr) frequencies (cm-1) -115.035417222 2.4*10-4 -36.6840, -0.0018, -0.0016, -0.0015, 15.6735, 28.1408, 343.6693 -115.035418323 2.2*10-6 -2.8167, -0.0017, -0.0014, -0.0006, 0.8876, 1.5931, 348.1832

This first conformer shows one larger imaginary frequency (-37 cm^-1), then three close-to-zero values, two small positive vibrations of 16 and 28 cm^-1 and finally one larger vibrational mode at 344 cm^-1. Is this a transition state or a local minimum? Even though we are clearly not having one large imaginary frequency in combination with sixe close-to-zero values as expected for a transition state, the six lowest vibrational frequencies are not all that close to zero. The second entry in the table above shows that a clear case can be made for a local minimum after reoptimization to a higher convergence criterion.

Geometry optimization of the eclipsed conformer of methanol using loose convergence criteria leads to a first structure whose vibrational frequency spectrum contains a number of imaginary vibrational frequencies (first entry in the following table):

E (Hartree) rms gradient (Hartree/Bohr) frequencies (cm-1) -115.033224934 9.8*10-4 -346.8098, -53.2629, -42.5572, -0.0020, -0.0011, -0.0011, 94.9730 -115.033252143 2.7*10-6 -351.2930, -2.1700, -2.0037, -1.4803, -0.0009, 0.0009, 0.0011

Reoptimization of the same structure with tighter convergence criteria leads to a much clearer result and we can recognize that the lowest true vibrational mode of the staggered conformer at 348 cm^-1 turns into an imaginary mode of 351 cm^-1 in the eclipsed structure. Both modes describe rotational motion around the central C-O bond.

In contrast to the analytical algorithm, failed numerical frequency calculations can be restarted with

freq=(Numerical,Restart)

Scaling Vibrational Frequencies
A final point concerns the accuracy obtained in vibrational frequency calculations as a function of the theoretical method used. In many cases, the experimentally measured values will differ significantly from those calculated theoretically. One source of discrepancy is that experimental values are often measured in solution while calculated values refer to the gas phase. A second source of error is the assumption of the harmonic approximation. There is, however, also a method-dependent difference between calculated and experimentally measured harmonic vibrational frequencies in the gas phase. For more economical methods such as HF or AM1, the deviations typically range around 10%. In order to account for this systematic deviation, all vibrational frequencies can be scaled (=multiplied with the same value) down to fit the experimental harmonic frequencies. Radom et al. have published a list of recommended scale factors for frequently used theoretical methods in J. Phys. Chem. 1996, 100, 16502. A small selection is:

method	scale factor	rms deviation (cm-1)
AM1	0.9532	126
PM3	0.9761	159
HF/6-31G(d)	0.8953	50
MP2(FC)/6-31G(d)	0.9427	61
QCISD(FC)/6-31G(d)	0.9537	37
BLYP/6-31G(d)	0.9945	45
B3LYP/6-31G(d)	0.9614	34

It is clear from these results that semiempirical methods such as AM1 and PM3 will, even after scaling, not be particularly reliable in predicting vibrational spectra. Rather good results can be obtained either from highly correlated, expensive methods such as QCISD, or from hybrid density functional calculations such as B3LYP.

Further Reading
A detailed discussion of further aspects of vibrational frequency calculations in Gaussian has been published by Joseph W. Ochterski (Gaussian Inc.) and can be found on the Gaussian home page: Vibrational Analysis in Gaussian.

last changes: 11.11.2004, HZ
questions & comments to: zipse@cup.uni-muenchen.de