Plotting theoretical infrared spectra
Published:
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 CO2 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 CO2 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
\[\mathcal{A}_\mathrm{e}= \int_\mathrm{band}{\alpha_\mathrm{e}\left(\widetilde{\nu}\right)\ \mathrm{d}\widetilde{\nu}} = \frac{1}{nl}\int_\mathrm{band}{A_\mathrm{e}\left(\widetilde{\nu}\right)\ \mathrm{d}\widetilde{\nu}} \\ = \frac{1}{nl}\int_\mathrm{band}{\ln\left[\frac{I_0(\widetilde{\nu})}{I(\widetilde{\nu})}\right]\left(\widetilde{\nu}\right)\ \mathrm{d}\widetilde{\nu}}\ \ \ \ \ \ \ \left[\mathrm{m \ mol^{-1}}\right]\]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
\[\mathcal{A}_e=\frac{N_Ad}{12\epsilon_0c^2}\left|\frac{\partial\mathbf{\mu}_\mathrm{elec.}}{\partial Q}\right|^2\cdot \mathrm{kg\left(u\right)\cdot \mathrm{D^2\left(C^2Å^2\right)} \ \ \ \ \ \ \ \ \ [m \ mol^{-1}] }\]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}}]$:
\[\mathrm{kg\left(u\right)} = \left(1.661\times{10}^{-27}\right)^{-1}\left[\mathrm{u\ kg}^{-1}\right]\\ \mathrm{D^2\left(C^2Å^2\right)}=\left(3.336\times10^{-20}\right)^2 \left[\mathrm{C^2Å^2D^{-2}}\right]\]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
\[L(x) = \frac{1}{\pi} \frac{\frac{1}{2}\Gamma}{\left(x-x_0\right)^2 + \left(\frac{1}{2}\Gamma\right)^2} \ \ \ \ [\mathrm{cm}] \\ \ \\ \int_{-\infty}^\infty L(x) \ \ \mathrm{d}x = 1\]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, its 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
\[\mathcal{A} \ \ [\mathrm{1000 \ cm \ mol^{-1}}] = \frac{100}{\ln{10}} \cdot \mathcal{A}_\mathrm{e} \ \ \ \ [\mathrm{km \ mol^{-1}}],\]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.