5.18. Vibrational Spectroscopy

The fundamental usage and several examples for computing vibrational frequencies in ORCA can be found here: Vibrational Frequencies. The following section is more spectroscopy focused.

5.18.1. IR Spectra

**There were significant changes in the IR printing after ORCA 4.2.1!**

IR spectral intensities are calculated automatically in frequency runs. Thus, there is nothing to control by the user. Consider the following job:

! OPT FREQ BP86 def2-SVP

*xyz 0 1
 O    0.000000    0.000000    0.611880
 C    0.000000    0.000000   -0.596849
 H    0.952616    0.000000   -1.209311
 H   -0.952616    0.000000   -1.209311
*

which gives the following output:

-----------
IR SPECTRUM
-----------

 Mode   freq       eps      Int      T**2         TX        TY        TZ
       cm**-1   L/(mol*cm) km/mol    a.u.
----------------------------------------------------------------------------
  6:   1146.68   0.000341    1.73  0.000093  (-0.000000 -0.009640  0.000000)
  7:   1224.67   0.002004   10.13  0.000511  ( 0.022596  0.000000  0.000000)
  8:   1485.77   0.001002    5.07  0.000211  ( 0.000000 -0.000000  0.014510)
  9:   1806.49   0.020286  102.51  0.003504  ( 0.000000 -0.000000  0.059197)
 10:   2769.13   0.014010   70.80  0.001579  ( 0.000000  0.000000  0.039734)
 11:   2812.52   0.039321  198.71  0.004363  ( 0.066052 -0.000000 -0.000000)

The first column (‘Mode’) labels vibrational modes that increase in frequency from top to bottom.” The next column provides vibrational frequencies. The molar absorption coefficient \(\varepsilon\) of each mode is listed in the “eps” column. This quantity is directly proportional to the intensity of a given fundamental in an IR spectrum, and thus it is used by the orca_mapspc utility program as the IR intensity.

The values under “Int” are the integrated absorption coefficient[1], and the “T**2” column lists the norm of the transition dipole derivatives, already including the vibrational part.

To obtain a plot of the spectrum, the orca_mapspc utility can be run calling the output file as:

orca_mapspc Test-Freq-H2CO.out ir -w25

or calling the Hessian file as:

orca_mapspc Test-Freq-H2CO.hess ir -w25

The basic options of orca_mapspc are listed below:

-w  : a value for the linewidth (gaussian shape, fwhm)
-x0 : start value of the spectrum in cm**-1
-x1 : end value of the spectrum in cm**-1
-n  : number of points to use

To see its options in detail, call orca_mapspc without any input. The above orca_mapspc runs of the H\(_2\)CO molecule provide Test-NumFreq-H2CO.out.ir.dat file that contains intensity and wavenumber columns. Therefore, this file can serve as input for any graph plotting program. The plot of the computed IR spectrum of the H\(_2\)CO molecule obtained with the above ORCA run is as given in Figure Fig. 5.54.

Fig. 5.54 The predicted IR spectrum of the H\(_2\)CO molecule plotted using the file generated by the orca_mapspc tool.

5.18.2. Overtones, Combination Bands and Near IR Spectra via NEARIR

Overtones and combination bands can also be incorporated to the computed IR or Near IR spectrum for completeness. The intensities of these bands are strongly dependent on anharmonic effects. ORCA can include these effects by means of the VPT2 approach [671]. The full cubic force field, anharmonic corrections to overtones and combination bands, and a broad range of methods are available in the orca_vpt2 module (see section Anharmonic Analysis and Vibrational Corrections using VPT2/GVPT2 and orca_vpt2).

In particular, the NEARIR keyword calls a simpler semidiagonal approach, including only two modes (\(i\) and \(j\), also refered as 2MR-QFF in [672, 673]) and force constants up to cubic order (\(k_{iij}\), \(k_{iji}\) and \(k_{iii}\)). For now, only the intensities are corrected for anharmonic effects - frequencies are not.

5.18.2.1. Overtones and Combination Bands

Since the calculation of these terms scale with \(N_{modes}^2\), it can quickly become too expensive, thus we use by default the semiempirical GFN2-xTB [318] to compute the energies and dipole moments necessary to the higher order derivatives (which can be changed later). To request this, simply add !NEARIR in the main input. An example input for computing the fundamentals of toluene using B2PLYP double-hybrid functional and for computing the anharmonics using XTB is as follows:

! TightOPT NumFREQ RI-B2PLYP def2-TZVP def2-TZVP/C RIJCOSX NEARIR
 
*xyzfile 0 1 toluene.xyz

Note

These anharmonic corrections are very sensitive to the geometry. Therefore, perform a conservative geometry optimization (at least TightOPT) whenever possible.

In the output, the characteristics of the regular IR spectrum are printed first. Then, the characteristics of overtones and combination bands are provided similarly to the fundamentals, as follows:

-------------------------------
OVERTONES AND COMBINATION BANDS
-------------------------------

Mode      freq       eps       Int      T**2        TX        TY        TZ
         cm**-1   L/(mol*cm)  km/mol    a.u.
------------------------------------------------------------------------------
6+6:      64.71    0.000994    5.02   0.004792  (-0.009428 -0.066232  0.017796)
6+7:      241.83   0.000022    0.11   0.000028  (-0.005268  0.000255  0.000638)
6+8:      375.36   0.000048    0.24   0.000040  (-0.000740  0.001917  0.006007)
6+9:      442.49   0.000000    0.00   0.000000  ( 0.000010  0.000001  0.000001)
6+10:     506.37   0.000003    0.01   0.000002  ( 0.001078 -0.000061  0.000799)
(...)

The “Mode” column shows the overtones, such as 6+6, and combination bands, such as 6+7 and 6+8. These new quantities are automatically detected and incorporated in the IR spectrum when the output file is called with the orca_mapspc utility as follows:

orca_mapspc toluene-nearir.out ir -w25

From the file orca_mapspc provided, the IR spectrum can be plotted as shown in Figure Fig. 5.55.

Fig. 5.55 Calculated and experimental infrared spectrum of toluene in gas phase. While the red plot includes only the fundamentals, the blue plot includes also overtones and combination bands. The grey dashed plot is the experimental gas-phase spectrum obtained from the NIST database. The theoretical frequencies were scaled following literature values [674]

“Benzene fingers”, i.e., overtones and combination bands of the ring, are recovered in the computed spectrum. Note that the frequencies were scaled using literature values [674], and are not yet corrected using VPT2.

5.18.2.2. Near IR Spectra

Let us simulate near IR spectrum of methanol in CCl\(_4\), as published by Bec and Huck [675], using B3LYP for fundamentals, XTB for overtones, and CPCM for solvation. The input is as follows:

! TightOPT FREQ B3LYP def2-TZVP RIJCOSX NEARIR CPCM(CCl4)

*xyz 0 1
 O      0.39517      4.38840     -0.00683
 C     -0.50818      3.29837      0.00221
 H     -0.11943      5.18771      0.19752
 H      0.03977      2.38083     -0.22470
 H     -1.27919      3.45664     -0.75583
 H     -0.96616      3.21170      0.99058
*

Calling the output with orca_mapspc by setting final point to about 8000 cm\(^{-1}\) in order to extend the spectrum to the near IR region, i.e.,

orca_mapspc toluene-nearir.out ir -w25 -x18000

one can simulate the spectrum from the generated “toluene-nearir.dat” file. As seen in Figure Fig. 5.56 the computed spectrum plotted with scaled computational frequencies (not yet corrected using VPT2) according to [674] agrees reasonably well with the experimental spectrum.

Fig. 5.56 Calculated and experimental near IR spectrum of methanol in CCl\(_4\). The blue plot is for overtones; the red plot is for combination bands; and the grey dashed plot is the experimental spectrum. Theoretical frequencies were scaled according to literature values [674].

5.18.2.3. Using Other Methods for the VPT2 Correction

To compute overtones with the method chosen for the calculation of the fundamentals, one needs only to set XTBVPT2 option in the %freq block to false, i.e.,

%freq XTBVPT2 False end

To set a different method for the calculation of overtones and combinations than used for the calculation of fundamentals, one needs first to perform a frequency calculation, then call the resulting Hessian file in %geom block, and activate the PRINTTHERMOCHEM flag (see section Thermochemistry for details), i.e.,

! BP86 def2-TZVP NEARIR CPCM(CCl4) PRINTTHERMOCHEM

%geom INHESSNAME "methanol.hess" end
%freq XTBVPT2 False end

*xyzfile 0 1 methanol_opt.xyz

In this example, the fundamental modes are read from the “methanol.hess” file, but the anharmonics and intensities of the overtones and combinations are computed using BP86. Any combination of methods, such as B3LYP/BP86 and B2PLYP/AM1, is allowed. Note that this description is an approximation to full VPT2 or GVPT2. For a more complete treatment, see the VPT2 module described in section Anharmonic Analysis and Vibrational Corrections using VPT2/GVPT2 and orca_vpt2.

By default, a step size of 0.5 in dimensionless normal mode unit is used during the numerical calculations. This can be changed by setting DELQ in the %freq block:

%freq
    XTBVPT2 False
    DELQ    0.1
end

The complete list of options related to VPT2 and in general frequency calculations can be found in Sec. Vibrational Frequencies.

5.18.3. Vibrational Circular Dichroism (VCD) Spectra

Vibrational circular dichroism spectrum calculations are implemented analytically at the SCF (HF or DFT) level following the derivation of Weigend and coworkers [676]. The basic usage is shown in the following example:

# AnFreq + doVCD triggers the VCD calculation
! AnFreq B3LYP def2-SVP

%freq 
  doVCD true 
end

*xyz 0 1
  C      1.231429   -0.226472   -0.084960
  C     -0.061893    0.507641    0.134338
  C     -1.358912   -0.147897    0.084831
  O     -0.902881    0.641038   -0.969176
  H      1.070541   -1.118875   -0.689778
  H      1.672013   -0.522768    0.869009
  H      1.946503    0.413187   -0.605194
  H      0.017832    1.411161    0.734623
  H     -1.417896   -1.212878   -0.118068
  H     -2.196737    0.255864    0.644375
*

Note that in addition to the Hessian, the VCD calculation requires the magnetic field response using GIAOs and the electric field response with the field origin placed at (0,0,0). The latter matches the hard-coded magnetic field gauge origin in the GIAO case and is necessary to ensure gauge-invariance of the results. ORCA does all of this automatically but it means that if VCD is requested together with electric and/or magnetic properties in the same job, the field origins cannot be changed.

Other keywords that influence the VCD calculation include GIAO_1el and GIAO_2el in %eprnmr and CutOffFreq in %freq. Note also that VCD cannot be computed with NumFreq.

5.18.4. Raman Spectra

Starting with ORCA 6.1, calculation of Raman intensities has been implemented using analytical geometric derivatives of the polarizability tensor and is now compatible with analytical frequency runs for SCF methods (HF, DFT). If a frequency run (!Freq) is combined with a polarizability calculation, the Raman spectrum will be automatically calculated.

Consider the following example:

! OPT RHF STO-3G TightSCF SmallPrint Freq

%elprop Polar 1 end

*xyz 0 1
 C    0.000000    0.000000   -0.533905
 O    0.000000    0.000000    0.682807
 H    0.000000    0.926563   -1.129511
 H    0.000000   -0.926563   -1.129511
*

The output provides the Raman scattering activity (in \({\text{Å}}^4\)/AMU) [677] and the Raman depolarization ratio of each mode:

--------------
RAMAN SPECTRUM
--------------

 Mode    freq (cm**-1)   Activity   Depolarization
-------------------------------------------------------------------
   6:      1278.43      0.007371      0.750000
   7:      1397.30      3.043852      0.750000
   8:      1767.07     16.391442      0.707250
   9:      2099.28      6.704492      0.075590
  10:      3498.60     38.681286      0.186510
  11:      3645.53     24.527417      0.750000

The ORCA run also generates a .hess file that includes polarizability derivatives and Raman activities. The effect of isotope substitution on the Raman activities can be computed using the .hess file.

As in the IR spectrum case, orca_mapspc provides a .dat file for plotting the computed Raman spectrum:

orca_mapspc Test-NumFreq-H2CO.out raman -w50

The Raman spectrum of H\(_2\)CO plotted by using the corresponding .dat file is as given in Figure Fig. 5.57.

Fig. 5.57 Calculated Raman spectrum of H\(_2\)CO at the STO-3G level plotted using the .dat generated by the orca_mapspc utility from numerical frequencies and Raman activities.

It is worth noting that Raman scattering activity \(S_i\) of each mode \(i\) is related to but not directly equal to the Raman intensity \(I_i\) of the corresponding mode, which is dependent on the excitation line \(\nu_0\) of the laser used in the Raman measurement(for Nd:YAG laser: \(\nu_0\) = 1064 nm = 9398.5 cm\(^{-1}\)). To obtain significantly better agreement between experimental and simulated Raman spectra, \(I_i\) of each mode needs to be computed with the following formula:

\[ I_i = \frac{f(\nu_0 - \nu_i)^4 S_i}{\nu_i [1 - \exp(-h c \nu_i / k T)]} \]

where \(f\) is a normalization constant common for all modes; \(h\), \(c\), \(k\), and \(T\) are Planck’s constant, speed of light, Boltzmann’s constant, and temperature, respectively.

Note

• It is still possible to calculate Raman intensities using numerical derivatives of the analytical polarizability, if the combination of !NumFreq and %elprop polar 1 end is specified. This is currently the only option for derivatives of the dynamical polarizability. In the static polarizability case, this is significantly less efficient and requires tighter accuracy settings than the analytical derivative.

• The analytical Raman intensities calculation is also implemented for implicit solvation, specified e.g. using !CPCM(solvent) in the simple input.

• The analytical Raman intensities make use of RI-J and/or COSX approximations in case these are used for the SCF. As a note of caution, we have observed that the error of the RI-J approximation may become larger with the default auxiliary basis in case you use a basis set with diffuse functions (e.g. def2-SVPD).

5.18.5. Resonance Raman Spectra

Resonance Raman spectra (NRVS) and excitation profiles can be predicted or fitted using the procedures described in section Simulation and Fit of Vibronic Structure in Electronic Spectra, Resonance Raman Excitation Profiles and Spectra with the orca_asa Program. An example for obtaining the necessary orca_asa input is described in section Absorption and Fluorescence Bandshapes using ORCA_ASA.

5.18.6. NRVS Spectra

The details of the theory and implementation of NRVS spectrum are as described in ref. [678, 679]. The NRVS spectrum of iron-containing molecules can be simply calculated calling .hess file of a previous frequency calculation with the orca_vib utility. The output file of this utility can then be called with orca_mapspc utility to produce a .dat file for plotting the spectrum:

orca_vib MyJob.hess > MyJob.vib.out
orca_mapspc MyJob.vib.out NRVS

For a the ferric-azide complex [678], the computed and experimental NRVS spectra are provided in Figure Fig. 5.58.

Fig. 5.58 Experimental (a, black curve), fitted (a, red) and simulated (b) NRVS spectrum of the Fe(III)-azide complex obtained at the BP86/TZVP level (T \(=\) 20 K). Bar graphs represent the corresponding intensities of the individual vibrational transitions. The blue curve represents the fitted spectrum with a background line removed.

As for the calculation of resonance Raman spectra described in section Simulation and Fit of Vibronic Structure in Electronic Spectra, Resonance Raman Excitation Profiles and Spectra with the orca_asa Program, the DFT estimations are usually excellent starting points for least-square refinements.

Below we describe the procedure for computing such NRVS spectra on the \({\text{Fe(SH)}}_{4}^{1-}\) complex with the BP86 functional, which typically provides good NRVS spectra. One needs first to optimize the geometry of the complex and compute its vibrational structure:

! OPT FREQ BP86 def2-TZVP TightSCF SmallPrint
 
*xyz -1 6
 Fe     -0.115452      0.019090     -0.059506
 S      -0.115452      1.781846      1.465006
 S      -0.115452     -1.743665      1.462801
 S      -1.908178     -0.072782     -1.518702
 S       1.560523      0.154286     -1.656664
 H       0.410700      2.760449      0.687716
 H      -0.674147     -2.708278      0.690223
 H      -2.905212      0.345589     -0.699907
 H       2.647892     -0.211681     -0.932926
*

Now run the orca_vib utility on the .hess file generated by this job to obtain an output file that can be used with orca_mapspc utility:

orca_vib Test-FeIIISH4-NumFreq.hess > Test-FeIIISH4-NumFreq.out
orca_mapspc Test-FeIIISH4-NumFreq.out NRVS

This orca_mapspc run generates Test-FeIIISH4-NumFreq.nrvs.dat file in the xy-format. This file contains phonon energy (\(x\), in cm\(^{-1}\)) and NRVS intensity (\(y\), in atomic units) and thus can be directly used for visualizing the spectrum.

The corresponding NRVS spectrum is given in Figure Fig. 5.59 together with the computational IR spectrum on the same frequency scale. NRVS reports the Doppler broadening of the Moessbauer signal due to resonant scattering of phonons (vibrations) dominated by the movements of Fe nuclei. This is a valuable addition to IR spectrum where the modes have very small intensities.

Fig. 5.59 (a) Theoretical IR spectrum of \({\text{Fe(SH)}}_{4}^{1-}\) and arrow-pictures of the highest intensity modes around the peak maxima. (b) The corresponding NRVS scattering spectrum.

5.18.7. Animation of Vibrational Modes

For describing how to animate vibrational modes and generate their “arrow-pictures”, let us perform a frequency calculation on H\(_2\)CO:

! OPT FREQ RHF STO-3G

*xyz 0 1
 C    0.000000    0.000000   -0.533905
 O    0.000000    0.000000    0.682807
 H    0.000000    0.926563   -1.129511
 H    0.000000   -0.926563   -1.129511
*

The output of this job provides vibrational characteristics:

Mode      freq       eps       Int    T**2         TX        TY        TZ
         cm**-1   L/(mol*cm) km/mol    a.u.
----------------------------------------------------------------------------
6:       1278.37   0.001222    6.18  0.000298  (-0.017272  0.000000  0.000000)
7:       1397.26   0.005844   29.53  0.001305  ( 0.000000  0.036128  0.000000)
8:       1767.02   0.000828    4.18  0.000146  (-0.000000  0.000000 -0.012089)
9:       2099.24   0.001668    8.43  0.000248  ( 0.000000 -0.000000  0.015749)
10:      3498.54   0.000356    1.80  0.000032  ( 0.000000 -0.000000 -0.005636)
11:      3645.47   0.003922   19.82  0.000336  (-0.000000  0.018322  0.000000)

This output can be directly opened with ChemCraft to visualize normal modes of H\(_2\)CO and to extract their arrow-pictures representing the direction of nuclear movements as shown in Figure Fig. 5.60. As an example, one can infer from this figure that the 1397 cm\(^{-1}\) mode corresponds to a rocking vibration.

Fig. 5.60 Normal modes of H\(_2\)CO with arrows indicating magnitude and direction of nuclear motions and the associated vibrational frequencies in cm\(^{-1}\) obtained from ORCA output file through the use of ChemCraft

In order to animate vibrational modes and to create their “arrow-pictures” by using free program packages like gOpenMol, the small utility program orca_pltvib can be used. This utility program generates a series of files from an ORCA output file of a frequency run, which can be openned with molecular visualization programs. The usage of orca_pltvib is as follows:

orca_pltvib Test-FREQ-H2CO.out [list of vibrations or all]

For example, let us animate the 1397 cm\(^{-1}\) mode labeled as 7:

orca_pltvib Test-FREQ-H2CO.out 7

This call will generate the Test-FREQ-H2CO.out.v007.xyz file. Open gOpenMol and read this file (Import->coords) in Xmol format. Then, go to the Trajectory->Main menu and import again the file in Xmol format. Now you are able to animate the mode. In order to generate a printable picture, press Dismiss and then type lulVectorDemo {4 0.1 black} into the gOpenMol command line window. The generated picture (see Figure Fig. 5.61) demonstrates that this mode corresponds to a rocking vibration.

Fig. 5.61 The 1397 cm\(^{-1}\) mode of the H\(_2\)CO molecule as obtained from the interface of ORCA to gOpenMol and the orca_pltvib tool to create the animation file.

The appearence of the arrows can also be modified as described in the web tutorial of gOpenMol.

5.18.8. Isotope Shifts

The calculated isotope shifts greatly aid in the identification of vibrations, the interpretation of experiments, and the assessment of the reliability of the calculated vibrational normal modes. It would be a very bad practice to recalculate the Hessian for investigating isotope shift since Hessian calculations are typically expensive, and the Hessian itself is independent of the masses. Below we describe how to find the isotope effect without recomputing the Hessian.

Let us suppose that you have calculated a Hessian as in the example discussed above, and you want to predict the effect of \(^{18}\)O substitution. In this case you can use the small utility program orca_vib. First of all you need to edit the masses given in the .hess file by hand. For the example given above, the .hess file is as follows:

$orca_hessian_file
......................
$hessian
12
... the cartesian Hessian in Eh/bohr**"

$vibrational_frequencies
12

...the vibrational frequencies (in cm-1) as in the output

$normal_modes
12 12
... the vibrational normal modes in Cartesian displacements
#
# The atoms: label  mass x  y  z
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#   Here we have changed 15.999 for oxygen into
#   18.0  in order to see the oxygen 18 effects
# !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
$atoms
4
 C     12.0110      0.000000     0.000000    -1.149571
 O     18.0000     -0.000000    -0.000000     1.149695
 H      1.0080     -0.000000     1.750696    -2.275041
 H      1.0080     -0.000000    -1.750696    -2.275041

$actual_temperature
  0.000000

$dipole_derivatives
12

... the dipole derivatives (Cartesian displacements)
#
# The IR spectrum
#  wavenumber T**2 TX TY  TY
#
$ir_spectrum
12
... the IR intensities

After changing the mass of O from 15.999 to 18.0 as shown above, let us call:

orca_vib Test-FREQ-H2CO.hess

This will recompute vibrational frequencies in the presence of \(^{18}\)O. Let us compare vibrational frequencies in the output of this run with the original frequencies in cm\(^{-1}\):

  Mode       H2C16O    H2CO18O    Shift
----------------------------------------- 
   6:        1284.36   1282.82    -1.54
   7:        1397.40   1391.74    -5.66
   8:        1766.60   1751.62   -14.98
   9:        2099.20   2061.49   -37.71
  10:        3499.11   3499.02    -0.09
  11:        3645.24   3645.24     0.00

Another way to analyze isotope shifts is to plot the two predicted spectra and then subtract one from the other. This will produce derivative-shaped peaks with zero crossings at the positions of the isotope-sensitive modes.

Note

In the presence of point charges and/or an external electric field, the translational and rotational symmetries of the system may be broken. In such cases, you may prefer NOT to project out the translational and rotational degrees of freedom of the Hessian. This can be achieved as:

orca_vib Test-FREQ-H2CO.hess -noproj

---

[1]

Explained in more detail by Neugbauer [677]