X

Search Results

Searching....

10.7 Harmonic Vibrational Analysis

10.7.1 Overview

(May 21, 2025)

Vibrational analysis is an extremely important tool for the quantum chemist, supplying a molecular fingerprint which is invaluable for aiding identification of molecular species in many experimental studies. Q-Chem includes a vibrational analysis package that can calculate vibrational frequencies and their infrared and Raman activities. 625 Johnson B. G., Florián J.
Chem. Phys. Lett.
(1995), 247, pp. 120.
Link
Vibrational frequencies are calculated by either using an analytic Hessian (if available) or else by numerical finite difference of the gradient or, as a last resort, via double finite-difference of single-point energies. (See Table 9.2 for the availability of analytic gradients and Hessians. The finite-difference step size is controlled by FDIFF_STEPSIZE.) The default setting in Q-Chem is to use the highest analytical derivative order available for the requested theoretical method. The performance of various ab initio theories in determining vibrational frequencies has been well documented. 930 Murray C. W. et al.
Chem. Phys. Lett.
(1992), 199, pp. 551.
Link
, 628 Johnson B. G., Gill P. M. W., Pople J. A.
J. Chem. Phys.
(1993), 98, pp. 5612.
Link
, 1170 Scott A. P., Radom L.
J. Phys. Chem.
(1996), 100, pp. 16502.
Link

When calculating analytic frequencies at the HF and DFT levels of theory, the coupled-perturbed SCF equations must be solved. This is the most time-consuming step in the calculation, and also consumes the most memory. The amount of memory required is 𝒪(N2M) where N is the number of basis functions, and M the number of atoms. This is an order more memory than is required for the SCF calculation, and is often the limiting consideration when treating larger systems analytically. Q-Chem incorporates a new approach to this problem that avoids this memory bottleneck by solving the CPSCF equations in segments. 689 Korambath P. P. et al.
Mol. Phys.
(2002), 100, pp. 1755.
Link
Instead of solving for all the perturbations at once, they are divided into several segments, and the CPSCF is applied for one segment at a time, resulting in a memory scaling of 𝒪(N2M/Nseg), where Nseg is the number of segments. This option is invoked automatically by the program.

After a vibrational analysis, Q-Chem computes useful statistical thermodynamic properties at standard temperature and pressure following the rigid-rotor-harmonic-oscillator (RRHO) approach. These include the zero-point vibration energy (ZPVE) and, translational, rotational and vibrational, entropies and enthalpies. Note: in the Q-Chem output the “total enthalpy” actually means the total enthalpy correction to the internal energy. One must add this “total enthalpy” to the internal energy to obtain the total enthalpy in common sense. In addition to these thermal corrections, Q-Chem also prints the interpolated vibrational entropy and enthalpy according to the quasi-RRHO (qRRHO) approach by Head-Gordon and co-workers, 794 Li Y.-P. et al.
J. Phys. Chem. C
(2015), 119, pp. 1840.
Link
which extends Grimme’s previous scheme 476 Grimme S.
Chem. Eur. J
(2012), 18, pp. 9955.
Link
to address low-frequency vibrations; Section 10.7.3.

10.7.1.1 Job Control

In order to carry out a frequency analysis users must at a minimum provide a molecule within the $molecule keyword and define an appropriate level of theory within the $rem keyword using the $rem variables EXCHANGE, CORRELATION (if required) (Chapter 4) and BASIS (Chapter 8). Since the default type of job (JOBTYPE) is a single point energy (SP) calculation, the JOBTYPE $rem variable must be set to FREQ.

It is very important to note that a vibrational frequency analysis must be performed at a stationary point on the potential surface that has been optimized at the same level of theory. Therefore a vibrational frequency analysis most naturally follows a geometry optimization in the same input deck, where the molecular geometry is obtained (see examples).

Users should also be aware that the quality of the quadrature grid used in DFT calculations is more important when calculating second derivatives. The default grid for some atoms has changed in Q-Chem 3.0 (see Section 5.5) and for this reason vibrational frequencies may vary slightly form previous versions. It is recommended that a grid larger than the default grid is used when performing frequency calculations.

The standard output from a frequency analysis includes the following.

  • Vibrational frequencies.

  • Raman and IR activities and intensities (requires $rem DORAMAN).

  • Atomic masses.

  • Zero-point vibrational energy.

  • Translational, rotational, and vibrational, entropies and enthalpies.

Several other $rem variables are available that control the vibrational frequency analysis. In detail, they are:

DORAMAN

DORAMAN
       Controls calculation of Raman intensities. Requires JOBTYPE to be set to FREQ
TYPE:
       LOGICAL
DEFAULT:
       FALSE
OPTIONS:
       FALSE Do not calculate Raman intensities. TRUE Do calculate Raman intensities.
RECOMMENDATION:
       None

VIBMAN_PRINT

VIBMAN_PRINT
       Controls level of extra print out for vibrational analysis.
TYPE:
       INTEGER
DEFAULT:
       1
OPTIONS:
       1 Standard full information print out. If VCI is TRUE, overtones and combination bands are also printed. 3 Level 1 plus vibrational frequencies in atomic units. 4 Level 3 plus mass-weighted Hessian matrix, projected mass-weighted Hessian matrix. 6 Level 4 plus vectors for translations and rotations projection matrix.
RECOMMENDATION:
       Use the default, unless the Hessian is desired.

CPSCF_NSEG

CPSCF_NSEG
       Controls the number of segments used to calculate the CPSCF equations.
TYPE:
       INTEGER
DEFAULT:
       0
OPTIONS:
       0 Determine the number of segments based on the memory request and MEM_TOTAL n User-defined. Use n segments when solving the CPSCF equations.
RECOMMENDATION:
       Use the default.

Example 10.30  EDF1/6-31+G* optimization followed by vibrational analysis. Doing the vibrational analysis at a stationary point is necessary for the results to be valid.

$molecule
   0  1
   O
   C  1  co
   F  2  fc  1  fco
   H  2  hc  1  hco  3  180.0

   co  =   1.2
   fc  =   1.4
   hc  =   1.0
   fco = 120.0
   hco = 120.0
$end

$rem
   JOBTYPE    opt
   METHOD     edf1
   BASIS      6-31+G*
$end

@@@

$molecule
   read
$end

$rem
   JOBTYPE    freq
   METHOD     edf1
   BASIS      6-31+G*
$end

10.7.1.2 Interpreting the Output

Numerical values in the following discussion correspond to Example 10.7.1.1 and the referenced partition functions come from the textbook by McQuarrie & Simon. Note that Q-Chem assumes T=298.15 K and P=1.00 atm by default; see Section 10.7.2 for instructions on how to modify these choices.

The quantity listed as zero-point vibrational energy (ZPVE),

Zero point vibrational energy: 12.695 kcal/mol

corresponds to ZPVE=12KhcνK. Within Q-Chem the vibrational entropy (Svib) is computed from the bottom of the well rather than from the zero-point level, and thus includes the ZPE. This corresponds to a vibrational partition function

qvib(T)=Ke-Θν,K/2T1-e-Θν,K/T (10.41)

where K indexes vibrational modes and ΘK=hνK/kB is the vibrational temperature of the Kth mode. This expression can be altered by a factor of ΘK/2 to start at the first vibrational level rather than the bottom of the well. The vibrational entropy is then

Svib=kBT(lnqvibT)V+kBlnqvib=kBTK[ΘK/Te-ΘK/T-1-ln(1-e-ΘK/T)], (10.42)

and for Example 10.7.1.1 it is reported as

Vibrational Entropy: 0.620 cal/mol.K

The vibrational enthalpy includes the ZPVE along with an additional temperature-dependent term:

Hvib=hcK[νK2+νKeΘK/T-1]. (10.43)

It is reported as

Vibrational Enthalpy: 12.839 kcal/mol

The rotational and translational enthalpy are multiples of RT,

Htotal=Hvib+Hrot+Htrans+RT, (10.44)

where the final contribution of RT comes from the definition H=U+PV. The quantity Htotal reported as

Total Enthalpy: 15.209 kcal/mol

Translational and rotational entropies can be computed from the corresponding ideal-gas partition functions.

For thermochemical calculations performed in implicit solvent (using models described in Section 11.2), there is a subtle point with some controversy attached. 564 Ho J., Klamt A., Coote M. L.
J. Phys. Chem. A
(2010), 114, pp. 13442.
Link
, 1116 Ribeiro R. F. et al.
J. Phys. Chem. B
(2011), 115, pp. 14556.
Link
, 205 Casasnovas R. et al.
Int. J. Quantum Chem.
(2014), 114, pp. 1350.
Link
, 548 Herbert J. M.
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021), 11, pp. e1519.
Link
This regards whether vibrational contributions to the free energy should be included or not. In nearly all cases, continuum solvation models are parameterized in order to reproduce experimental free energies of solvation (ΔsolvG) using rigid gas-phase geometries for the solute molecules. There is a potential double-counting problem if -TSvib from a vibrational frequency calculation in implicit solvent is added, which might be avoided by instead using gas-phase harmonic frequencies for the solute. 564 Ho J., Klamt A., Coote M. L.
J. Phys. Chem. A
(2010), 114, pp. 13442.
Link
However, the difference between these procedures is negligible (0.2 kcal/mol). for the small-molecule data sets that are used to train implicit solvent models, 1116 Ribeiro R. F. et al.
J. Phys. Chem. B
(2011), 115, pp. 14556.
Link
and only by using solution-phase vibrational frequencies can one obtain corrections to Svib arising from solvation-induced changes in geometry, which might be significant for larger molecules. Note that Q-Chem’s harmonic frequency engine assumes a gas-phase molecule (even in the presence of continuum solvent), so that rotational and translational contributions to the free energy are printed out in every case. These should not be included in solution-phase free energy differences.

10.7.1.3 Finite-Difference Frequency Calculation with Molecular Segments

For electronic structure methods whose analytic 2nd nuclear derivatives are not implemented in Q-Chem, finite-difference hessian calculations based on nuclear gradients will be performed. The Hessian matrix elements are calculated as

Hij=2Exixj=gi(xj+Δxj)-gi(xj-Δxj)2Δxj (10.45)

where gi=E/xi. In total, 6Natom gradient calculations at displaced structures are needed to construct the Hessian.

By default, the 6Natom gradient calculations are performed in serial, which is extremely inefficient for molecules of a large number of atoms. Taking advantage of the fact that these finite-difference calculations are independent of each other, Q-Chem provides the option to perform these calculations in molecular segments (i.e., atom blocks). The user will need to specify the number of atoms in each segment using $rem variable FD2ND_BLOCK_SIZE. The total number of segments will then be Natom/blocksize+1 (note: the number of atoms in the last segment is allowed to be smaller than the user-specified block size). The user will then need to set up and execute a series of input with FD2ND_BLOCK_INDEX set to differnt values (starting from 0). Note that all these jobs should share the same value of FD2ND_BLOCK_SIZE. On a large computer cluster, these jobs can easily be executed simultaneously. These jobs will produce text files grdntDisp.0, grdntDisp.1, … and dipoleDisp.0, dipoleDisp.1… under the working directory. Note that IDERIV = 1 needs to be added if it is not set automatically.

After all the calculations for molecule segments finish, the user will need to perform an additional job with FD2ND_BLOCK_INDEX set to -1. This job will collect the results for all the segments (from the text files in the working directory), construct the Hessian, and complete the vibrational frequency analysis. This approach will greatly speed up the finite-difference frequency calculations for large molecules.

FD2ND_BLOCK_SIZE

FD2ND_BLOCK_SIZE
       Controls the number of atoms in each segment of finite-difference Hessian calculation.
TYPE:
       INTEGER
DEFAULT:
       Uninitialized
OPTIONS:
       n Having n atoms in each segment
RECOMMENDATION:
       The value can be estimated based on the desired number of segments. It should be significantly smaller than the total number of atoms in the molecule.

FD2ND_BLOCK_INDEX

FD2ND_BLOCK_INDEX
       The flag that species which atom block is being processed in a given finite-difference calculation.
TYPE:
       INTEGER
DEFAULT:
       Uninitialized
OPTIONS:
       0,1,,Nfrag-1 The index of the first fragment is 0 and the last should be Nfrag-1 -1 Indicator of the job that wraps up the result.
RECOMMENDATION:
       None

Example 10.31  Finite-difference frequency calculation for a water dimer with two segments. In a production calculation, the first two jobs should be completed first as two separate calculations, and the 3rd job will be executed afterwards as a wrap-up.

$molecule
0 1
O    -1.48239   -0.12728    0.00000
H    -1.90743    0.74236    0.00000
H    -0.52048    0.04679    0.00000
O     1.37121    0.12361    0.00000
H     1.71353   -0.35106    0.77313
H     1.71353   -0.35106   -0.77313
$end

$rem
JOBTYPE  FREQ
METHOD  B3LYP
BASIS  6-31+G(D)
SYM_IGNORE  TRUE
THRESH      14
SCF_CONVERGENCE  8
IDERIV      1
FD2ND_BLOCK_SIZE  3
FD2ND_BLOCK_INDEX 0
$end

@@@

$molecule
read
$end

$rem
JOBTYPE  FREQ
METHOD  B3LYP
BASIS  6-31+G(D)
SYM_IGNORE  TRUE
THRESH      14
SCF_CONVERGENCE  8
IDERIV      1
FD2ND_BLOCK_SIZE  3
FD2ND_BLOCK_INDEX 1
$end

@@@

$molecule
read
$end

$rem
JOBTYPE  FREQ
METHOD  B3LYP
BASIS  6-31+G(D)
SYM_IGNORE  TRUE
THRESH      14
SCF_CONVERGENCE  8
IDERIV      1
FD2ND_BLOCK_SIZE  3
FD2ND_BLOCK_INDEX -1
$end