(sec:spectroscopyproperties.vib)= # Vibrational Spectroscopy The fundamental usage and several examples for computing vibrational frequencies in ORCA can be found here: {ref}`sec:structurereactivity.frequencies`. The following section is more spectroscopy focused. ```{index} IR, Vibrational Spectroscopy, Infrared Spectra ``` (sec:spectroscopyproperties.vib.ir)= ## 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: ```orcainput ! 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: ```orca ----------- 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 orca_mapspc Test-Freq-H2CO.out ir -w25 ``` or calling the Hessian file as: ```orca orca_mapspc Test-Freq-H2CO.hess ir -w25 ``` The basic options of `orca_mapspc` are listed below: ```orca -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 {numref}`fig:637`. (fig:637)= :::{figure} ../../images/637.* The predicted IR spectrum of the H$_2$CO molecule plotted using the file generated by the `orca_mapspc` tool. ::: (sec:spectroscopyproperties.vib.nearir)= ## Overtones, Combination Bands and Near IR Spectra via NEARIR ```{index} Overtones, Combination Bands, 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 {cite}`Bloino2014`. 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 {ref}`sec:spectroscopyproperties.vpt2`). In particular, the `NEARIR` keyword calls a simpler semidiagonal approach, including only two modes ($i$ and $j$, also refered as 2MR-QFF in {cite}`yagi_ab_2004,barnes_fast_2016`) 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.** ### 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 {cite}`grimme2017xtb` 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: ```orcainput ! 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: ```orca ------------------------------- 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 orca_mapspc toluene-nearir.out ir -w25 ``` From the file `orca_mapspc` provided, the IR spectrum can be plotted as shown in Figure {numref}`fig:nearir_toluene`. (fig:nearir_toluene)= :::{figure} ../../images/nearir_toluene.* 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 {cite}`keshfreq2015` ::: "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 {cite}`keshfreq2015`, and are not yet corrected using VPT2. ```{index} Near IR (NIR) ``` ### Near IR Spectra Let us simulate near IR spectrum of methanol in CCl$_4$, as published by Bec and Huck {cite}`bec_breakthrough_2019`, using B3LYP for fundamentals, XTB for overtones, and CPCM for solvation. The input is as follows: ```orcainput ! 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 orca_mapspc toluene-nearir.out ir -w25 -x18000 ``` one can simulate the spectrum from the generated "toluene-nearir.dat" file. As seen in Figure {numref}`fig:nearir_methanol` the computed spectrum plotted with scaled computational frequencies (not yet corrected using VPT2) according to {cite}`keshfreq2015` agrees reasonably well with the experimental spectrum. (fig:nearir_methanol)= ```{figure} ../../images/nearir_methanol.* 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 {cite}`keshfreq2015`. ``` ### 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., ```orca %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 {ref}`sec:structurereactivity.thermochemistry` for details), i.e., ```orca ! 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 {ref}`sec:spectroscopyproperties.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: ```orca %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. {ref}`sec:structurereactivity.frequencies`. (sec:spectroscopyproperties.vcd.typical)= ## Vibrational Circular Dichroism (VCD) Spectra ```{index} Vibrational Circular Dichroism ``` Vibrational circular dichroism spectrum calculations are implemented analytically at the SCF (HF or DFT) level following the derivation of Weigend and coworkers {cite}`Reiter2017`. The basic usage is shown in the following example: ```{literalinclude} ../../orca_working_input/VCD.inp :language: orca ``` 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`. (sec:spectroscopyproperties.vib.raman)= ## Raman Spectra ```{index} 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: ```orcainput ! 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) {cite}`Neugebauer2002` and the Raman depolarization ratio of each mode: ```orca -------------- 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 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 {numref}`fig:638`. (fig:638)= :::{figure} ../../images/638.* 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`). ::: (sec:spectroscopyproperties.vib.resonanceraman)= ## Resonance Raman Spectra ```{index} Resonance Raman Spectra ``` Resonance Raman spectra (NRVS) and excitation profiles can be predicted or fitted using the procedures described in section {ref}`sec:spectroscopyproperties.asa`. An example for obtaining the necessary orca_asa input is described in section {ref}`sec:spectroscopyproperties.asa.bandshapes`. (sec:spectroscopyproperties.vib.nrvs)= ## NRVS Spectra ```{index} NRVS Spectra ``` The details of the theory and implementation of NRVS spectrum are as described in ref. {cite}`petrenko2007hyperfine,debeer2007vibrational`. 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 orca_vib MyJob.hess > MyJob.vib.out orca_mapspc MyJob.vib.out NRVS ``` For a the ferric-azide complex {cite}`petrenko2007hyperfine`, the computed and experimental NRVS spectra are provided in Figure {numref}`fig:639`. (fig:639)= :::{figure} ../../images/639.* 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 {ref}`sec:spectroscopyproperties.asa`, 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: ```orcainput ! 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 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 {numref}`fig:639_1` 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:639_1)= :::{figure} ../../images/639_1.* (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. ::: (sec:spectroscopyproperties.vib.animation)= ## 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: ```{literalinclude} ../../orca_working_input/C05S13_213.inp :language: orca ``` The output of this job provides vibrational characteristics: ```orca 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 {numref}`fig:VibModes`. As an example, one can infer from this figure that the 1397 cm$^{-1}$ mode corresponds to a rocking vibration. (fig:VibModes)= :::{figure} ../../images/VibModes.* 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 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 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 {numref}`fig:640`) demonstrates that this mode corresponds to a rocking vibration. (fig:640)= :::{figure} ../../images/640.* 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`. (sec:spectroscopyproperties.vib.isotopeshifts)= ## Isotope Shifts ```{index} Isotope Shifts in Vibrational Spectra ``` 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 $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 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}$: ```orca 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 orca_vib Test-FREQ-H2CO.hess -noproj ``` ::: [^1]: Explained in more detail by Neugbauer {cite}`Neugebauer2002`