Plotting theoretical infrared spectra

I recently had to plot some theoretical infrared spectra using the output of a DFT optimisation and frequency calculation performed in Gaussian and didn’t want to have to plot everything in Gaussview. Here’s how to plot exactly the same spectrum as Gaussview, and hopefully a clear explanation of what Gaussian means by IR Intensities.

Running the calculation

Here’s an example optimisation and frequency (opt+freq) calculation for CO₂ using DFT (B3LYP XC Functional), and a reasonable basis set.

$RunGauss 
#p opt freq=HPModes b3lyp int=grid=ultrafine

DFT B3LYP OPT FREQ This-is-a-comment-line

0 1
O1     0      0.0000000   0.0000000   1.1600000
C1     0      0.0000000   0.0000000   0.0000000
O2     0      0.0000000   0.0000000  -1.1600000

C 0 
cc-pVTZ
**** 
O 0 
cc-pVTZ

I’m not going to explain how to run a calculation here in detail, but you’ll need freq in your routecard as above, and I would recommend freq=HPModes to print vibrational mode information with a reasonable number of significant figures.

Understanding the result

Gaussian will print the information you need to the output (.log) file in a section with the following header

 Harmonic frequencies (cm**-1), IR intensities (KM/Mole), Raman scattering
 activities (A**4/AMU), depolarization ratios for plane and unpolarized
 incident light, reduced masses (AMU), force constants (mDyne/A),
 and normal coordinates:

When using freq=HPModes as above, it will print this header twice, where the first section contains more precise displacements.

Sticking with the CO₂ example, the .log file contains

                          PIU       PIU       SGG       SGU
       Frequencies ---   487.4740  487.4740 1269.3077 2381.0728
    Reduced masses ---    12.8774   12.8774   15.9949   12.8774
   Force constants ---     1.8029    1.8029   15.1833   43.0153
    IR Intensities ---     9.2871    9.2871    0.0000   88.9346
 Coord Atom Element:
   1     1     8          0.00041  -0.33138  -0.00000  -0.00000
   2     1     8         -0.33138  -0.00041   0.00000  -0.00000
   3     1     8          0.00000  -0.00000   0.70711  -0.33138
   1     2     6         -0.00111   0.88339   0.00000   0.00000
   2     2     6          0.88339   0.00111  -0.00000   0.00000
   3     2     6         -0.00000  -0.00000  -0.00000   0.88339
   1     3     8          0.00041  -0.33138  -0.00000  -0.00000
   2     3     8         -0.33138  -0.00041   0.00000  -0.00000
   3     3     8         -0.00000   0.00000  -0.70711  -0.33138

We only need the Frequencies and IR Intensities which are in units of cm^-1 and km mol^-1 respectively.

       Frequencies ---   487.4740  487.4740 1269.3077 2381.0728
    IR Intensities ---     9.2871    9.2871    0.0000   88.9346

The former are easy to understand, they’re the fundamental frequency of each vibrational mode computed within the harmonic approximation. The latter are more complicated, and are the integrated napierian (natural logarithmic) absorbances of each mode.

Experimentally, this is the integral of an absorption band

Where $n$ is concentration $[\mathrm{mol \ m^{-3}]}$, $l$ is path length $[\mathrm{m}]$, $A_\mathrm{e}$ is the napierian absorbance, $\alpha_\mathrm{e}$ is the napierian absorption coefficient $[\mathrm{m^2 \ mol^{-1}}]$, $I$ and $I_0$ are the transmitted and incident radiation intensity respectively, and $\widetilde{\nu}$ is the wavenumber $[\mathrm{cm}^{-1}]$.

This quantity is calculated theoretically as

Where $N\mathrm{_A}$ is Avogadro’s number $\mathrm{[mol^{-1}]}$, $d$ is the degeneracy of the mode (assumed to be 1 for all modes), $c$ is the speed of light $\mathrm{[m \ s^{-1}]}$, and $\epsilon_0$ is the permittivity of free space $\mathrm{[F \ m^{-1}]}$. The derivative is that of the electric dipole moment with respect to displacement along the $x$, $y$, and $z$ components of the mode vector $\mathrm{[D \ Å^{-1} u^{-1/2}]}$. The final two terms explicitly detail the conversion factors required to arrive at the final units of $[\mathrm{m \ mol^{-1}}]$:

So when Gaussian prints the IR Intensities of the modes, it means $\mathcal{A}_\mathrm{e}$.

Plotting a spectrum

Gaussview uses an area-normalised Lorentzian lineshape $L(x)$ for each vibrational mode

where $\Gamma$ is the full-width-at-half-maximum (FWHM) of the curve, $x$ is the wavenumber of the incident radiation, and $x_0$ is the wavenumber of the vibrational mode.

The area of this function is then scaled to match the linear integrated absorption coefficent $\mathcal{A}$.

To do this, first divide Gaussian’s value of $\mathcal{A}_\mathrm{e} \ \ [\mathrm{km \ mol^{-1}}]$ by $\ln(10)$ to move from logarithmic to linear absorption coefficient. Then, it’s useful (see below) to multiply by 100 to change from $[\mathrm{km \ mol^{-1}}]$ to $[\mathrm{1000 \ cm \ mol^{-1}}]$. This may seem confusing, but if you realise that $[\mathrm{km \ mol^{-1}}] = [\mathrm{1000 \ m \ mol^{-1}}]$ then it might make this step clearer. To sum up

The multiplication of $L(x) \ [\mathrm{cm}]$ by $\mathcal{A} \ \ [\mathrm{1000 \ cm \ mol^{-1}}]$ gives a function with units of $[\mathrm{1000 \ cm^{2} \ mol^{-1}}] = [\mathrm{L \ mol^{-1} \ cm^{-1}}] = [\mathrm{M^{-1} \ cm^{-1}}]$, i.e. common units of the linear molar extinction coefficient $\epsilon$.

To combine all of the modes as a single spectrum, simply sum $L(x) \cdot \mathcal{A}$ for each mode, and plot against wavenumber.

Addendum - understanding Dipole Derivatives from Gaussian

If you want to obtain dipole derivatives with respect to normal modes using iop(7/33=1), then Gaussian will print these in units of $\mathrm{\sqrt{km\,mol^{−1}}}$ (without saying so). Various forum posts and online resources tell you to divide by the square root of 42.2561 to get a dipole derivative in more reasonable units (i.e. in the above definition), but they often don’t tell you where this number comes from.

The Gaussian website also lists the conversion factor

1 Debye2-angstrom-2-amu-1 (IR intensity unit) = 42.2561 km-mol-1

This value comes from evaluating the constant terms in the definition of $\mathcal{A}_\mathrm{e}$.

All of this goes to say, be careful when using the output of Gaussian!

Addendum - Dipole Derivatives in Cartesian basis

You can also obtain dipole derivatives with respect to normal modes without iop(7/33=1) by looking in the .fchk file for your calculation (obtained with Gaussian’s formchk utility). You will quickly notice, however, that these numbers do not match those in the output file, as they are the derivatives with respect to single atom displacements along $x$, $y$, and $z$.

To obtain the derivatives with respect to normal modes, you first must generate the matrix of cartesian displacements $\mathbf{L_\mathrm{CART}}$, as detailed in the Gaussian vibrational analysis whitepaper. You can take this matrix from the output file but beware, these are scaled by the square root of the reduced mass and so this scaling must be undone to get the true $\mathbf{L_\mathrm{CART}}$.

Then, for each mode, multiply each element of $\mathbf{L_\mathrm{CART}}$ by the per-atom dipole derivative vector from the .fchk file, and sum up all contributions in each direction to give a single three element vector describing the dipole derivative with respect to a given mode in $x$, $y$, and $z$. Finally, convert these values from $\left[\mathrm{e\,u^{-\frac{1}{2}}}\right]$ to $\left[\sqrt{\mathrm{km\, mol^{-1}}}\right]$ by multiplying by 31.2231, this then gives the same values as those printed by the IOp command.