IR spectra, (resonance) Raman, VROA, VCD

In ADF infrared and Raman spectroscopy can be studied for molecular vibrations. In the Born-Oppenheimer and harmonic approximations the vibrational frequencies are determined by the normal modes corresponding to the molecular electronic ground state potential energy surface. In resonance Raman spectroscopy the molecule is excited to near one of its electronic excited states, to improve the sensitivity compared to traditional Raman spectroscopy. Vibrational circular dichroism (VCD) is the differential absorption of left and right circularly polarized infrared light by vibrating molecules.

IR spectra

The IR frequencies can be calculated with the FREQUENCIES sub-block of the key GEOMETRY (numerical frequencies),

GEOMETRY
  FREQUENCIES
  End
END

or with the block key ANALYICALFREQ (analytical frequencies),

ANALYICALFREQ
END

These keys are described more extensively here: IR Frequencies.

Raman scattering

Raman scattering intensities and depolarization ratios for all molecular vibrations at a certain laser frequency can be calculated in a single run. The run type must be Frequencies, which is arranged with the FREQUENCIES subkey of the key GEOMETRY (numerical frequencies), or with the block key ANALYICALFREQ (analytical frequencies), see IR Frequencies.

Here the RESPONSE key is used to specify that Raman intensities are computed. In case of spin-orbit coupling or hybrids one needs to use the AORESPONSE key. The frequency dependent Raman scattering can be calculated for 1 laser frequency at the time.

RESPONSE
  RAMAN
  Nfreq 1
  FrqBeg Laserfreq
END
Frequencies or wavelengths
The number of frequencies Nfreq should be 1. With subkey Frqbeg the value of the Laser frequency value (Laserfreq) can be given. The unit of Frqbeg is eV.

For static Raman scattering (\(\omega = 0\)) use:

RESPONSE
 RAMAN
END

The Raman scattering calculation is very similar to an IR intensity calculation. In fact, all IR output is automatically generated as well. At all distorted geometries the dipole polarizability tensor is calculated. This is very time-consuming and is only feasible for small molecules. More details on the RESPONSE key can be found here.

There are a few caveats:

  • Numerical integration (BeckeGrid) accuracy must be high
  • A calculation in which only a subset of the atoms is displaced is not possible for Raman calculations.
  • For good results, a well converged (with the same basis and functional) equilibrium geometry must be used.

Because of this last point, it is wise to always start the RAMAN calculation with a TAPE13 restart file from a previous geometry optimization with the same basis, accuracy parameters, and density functional.

Atomic coordinate displacements in a RAMAN calculation must be Cartesian, not Z-matrix. Furthermore, the current implementation does not yet support constrained displacements, i.e. you must use all atomic coordinate displacements. However, one can calculate Raman for selected frequencies, see next section.

The alternative Raman implementation with the AORESPONSE offers some unique features like lifetime options, hybrids support, spin-orbit coupling.

AORESPONSE
 RAMAN
END

Raman Intensities for Selected Frequencies

The RAMANRANGE keyword can be used to calculate Raman intensities for a range of frequencies only. Recommended to be used in the case one use a t21 as a restart file, which has frequencies on them. Using this option is a fast alternative for the existing method of calculating Raman intensities. By default RAMANRANGE uses the RESPONSE key, however, in case of hybrids or spin-orbit coupling RAMANRANGE will use the AORESPONSE key. Raman icw hybrids icw spin-orbit coupling is not supported in ADF. Note: while RAMANRANGE works in combination with SYMMETRY NOSYM, it might not if the user explicitly specifies any other symmetry via the SYMMETRY key. The input keyword is as follows:

RAMANRANGE low high {NUM=num DISRAD=disrad}
low, high
Two values defining an interval of frequencies to calculate the Raman intensities for. The Raman intensities are calculated by numerical differentiation of the polarizability tensor. Only frequencies within the interval that are known to be Raman-active will be included. This means that 2**N* single-point TDDFT calculations will be performed, where N is the number of Raman-active frequencies within the range. Imaginary frequencies are specified using negative values, which is consistent with the notation adopted within ADF.
num
Num is an integer number specifying how many points are to be used for numerical differentiation: 2, 4, or 6. The default value is 2.
disrad
Disrad specifies the step size (in Angstrom) to be made in each direction for 2-point differentiation. For the 4-point differentiation the maximum deviation from the equilibrium geometry will be twice as large as for the 2-point one and so on. The default value for disrad is the same as used for numerical frequencies.

Main advantages of the method:

  • Full symmetry at each displaced geometry is used, which does not only speeds the calculation up but also makes it more accurate.
  • Only Raman-active modes are included in the calculation, which may save a lot of time for molecules with symmetry.
  • There is no need to recalculate frequencies if you already have a t21 file with them as it can be used as a restart file.

For static Raman scattering (\(\omega\) = 0) one does not need to add the RESPONSE (or AORESPONSE) block key. However, for the calculation of the frequency dependent Raman scattering the following RESPONSE block key is needed in the input:

RESPONSE
 RAMAN
 Nfreq 1
 FrqBeg Laserfreq
END
Frequencies or wavelengths
The number of frequencies Nfreq should be 1. With subkey Frqbeg the value of the Laser frequency value (Laserfreq) can be given. The unit of Frqbeg is eV.

For the AORESPONSE key (see also next section), this would be:

AORESPONSe
 RAMAN
 Frequency 1 freq1 units
 {LIFETIME width}
END
Frequency 1 freq1 units
To calculate time-dependent properties, one needs to specify the frequency of perturbation field. Here 1 frequency freq1 is used. The last item on the line specifies the units and is one of EV, HARTREE, ANGSTROM.
LIFETIME width`
Optional. See also next section, and the AORESPONSE documentation.

Resonance Raman: excited-state finite lifetime

In this method (Ref. [1]) the resonance Raman-scattering (RRS) spectra is calculated from the geometrical derivatives of the frequency-dependent polarizability. The polarizability derivatives are calculated from resonance polarizabilities by including a finite lifetime (phenomenological parameter) of the electronic excited states using time-dependent density-functional theory.

It is similar to the simple excited-state gradient approximation method (see next section) if only one electronic excited state is important, however, it is not restricted to only one electronic excited state. In the limit that there is only one possible state in resonance the two methods should give more or less the same results. However, for many states and high-energy states and to get resonance Raman profiles (i.e., Raman intensities as a function of the energy of the incident light beam) this approach might be more suitable. The resonance Raman profiles in this approach are averaged profiles since vibronic coupling effects are not accounted for. At the moment this method needs numerically calculated frequencies in Cartesian coordinates. The method described in the next section can use a basis of normal coordinates rather than Cartesian coordinates, so that in that method the calculation can be restricted to a couple of modes.

GEOMETRY
 FREQUENCIES
 End
END
AORESPONSE
 RAMAN
 FREQUENCY  1 freq1 units
 LIFETIME width
END

The RRS is not calculated if one uses symmetric displacements in the numerically calculated frequencies or if one uses analytically calculated frequencies. However, one can use RAMANRANGE, see previous section. The documentation of the AORESPONSE key explains in more detail the meaning of the subkeywords in the block key AORESPONSE, which are required to calculate RRS. Similarly to the RESPONSE-Raman module, the AOREPONSE-Raman only works with one frequency.

Resonance Raman: excited-state gradient

According to a the time-dependent picture of resonance-Raman (RR) scattering the relative intensities of RR scattering cross sections are, under certain assumptions, proportional to the square of the excited-state energy gradients projected onto the ground-state normal modes of the molecule (see Ref. [2]). For an alternative implementation of RR scattering using a finite lifetime of the excited states, and a discussion of some of the differences, see the previous section.

For ADF two implementations exist that use the excited state gradient, also called vertical gradient Franck-Condon (VG-FC) Resonance Raman:

In practice both methods need the ground-state normal modes as input. The AMS implementation of RR also needs as input the excited state gradient at the ground state geometry. In the Vibron implementation of RR the excited state gradient is calculated using finite displacements, which can be relatively expensive. The AMS implementation of RR can also be used via the ADF-GUI, and is relatively easy to use.

Note

For VG-FC Resonance Raman we recommend using the AMS VG-FC Resonance Raman implementation with ADF used as an AMS engine.

VIBRON module

The excited-state gradients which are needed in this method can be computed numerically by ADF’s VIBRON module, which is invoked by selecting the VIBRON runtype in the GEOMETRY block key, the use of the VIBRON block key and the EXCITATION block key:

GEOMETRY
  VIBRON
END
VIBRON
 NMTAPE filename
 RESRAMAN
 {...
 ..}
END
EXCITATIONS
  LOWEST nlowest
END

The VIBRON module always requires an EXCITATIONS input block, in which the total number of excited states to be calculated must be specified. NMTAPE is the only obligatory keyword for the VIBRON module. It specifies the name of a TAPE21 file from a previous frequency calculation. This TAPE21 file is needed to read the normal modes w.r.t. which the derivatives are computed. I.e., a separate frequency calculation must be carried out first. The second subkeyword RESRAMAN invokes the resonance Raman calculation. Note that the VIBRON module is not suited for open-shell TDDFT.

Resonance Raman for several excited states

The numerical evaluation of resonance Raman intensities has the advantages that

  • Relative Intensities can be computed for several excited states at a time, since all excitation energies are determined simultaneously.
  • Intensities can be computed for a selected no. of modes.

The intensities calculated for two different states cannot directly be compared, since the excited-state gradients only provide relative intensities for each excited state in resonance. For RR intensities from different excited states, also other quantities play a role. The most important one is the transition dipole moment to the excited state in resonance, which enters the intensity expression with to the fourth power.

Restrictions: (avoided) crossings between excited-states

The numerical calculation of excited-state gradients has a number of advantages, but also a possible problem: If the step size is chosen too large, or if there are close-lying excited-states, then the order of the excited states can change. For such cases, the excited-state gradient method to estimate relative RR intensities is not reliable: If states with different electronic character (but of the same symmetry) are close in energy, this will cause an avoided crossing. If the numerical derivatives are, in this case, computed w.r.t. the adiabatic states, they will probably not reflect the true situation. Especially if the coupling matrix elements between the two excited states is small, the spectroscopic properties often behave as if there is no avoided crossing, i.e., according to the diabatic states. Such cases should be handled with extreme care, since it is often not possible in advance to see whether the adiabatic or the diabatic picture should be invoked.

Because of the possible (avoided) crossings the user must make sure that always enough excited states are calculated to include the state(s) of interest. E.g., if the resonance Raman intensities are required for the first excited state, also some higher excited states have to be included in the excitation calculation, as the first excited state at the ground-state equilibrium might be higher in energy for displaced structures.

Restrictions: results not trustworthy for higher excited states

Users should be aware of another technical point: Excited states are usually calculated from a Davidson diagonalization procedure, i.e., only a small number of eigenvalues and eigenvectors describing the lowest excitations are obtained. During finite displacements, some of the higher calculated states might leave the calculated energy window, while others enter it. Hence, the character of some of the higher calculated states can change. In such a case, the numerical differentiation based on a (simple) diabatic pictures will fail for the higher states, since no mapping between the excited states for reference (equilibrium) and displaced structure can be carried out.

The solution is rather simple: Users should always ask for more excited states than they are actually interested in, and discard the data for higher states, in particular for those which could not successfully be mapped for displaced structures (look for messages in the output like ‘State No. X cannot be expressed in terms of reference states’).

For advanced users it should be mentioned that it is possible to set an energy window within the range of states calculated, and only the states within this energy window will be taken into account in the evaluation. See the subkey ELTHRESH and EUTHRESH of the block key VIBRON.

Furthermore, it is possible to pick out certain states from this energy window, and only perform the mapping (and diabatization, if requested) and differentiation for them. See the subkey SELSTATE of the block key VIBRON.

Advanced Restarts

In some cases, it only becomes obvious which states have to be included in a (simple) diabatization after excitation energies for all displaced structures are calculated. Therefore, the selected states and the energy window settings can also be adjusted in a restart (with the usual restart key) after all single-point calculations are done. However, this is only possible if all raw data are saved to TAPE21, which might be an enormous amount of data.

Therefore, the user has to specify the subkeyword SAVERAWDAT of the key VIBRON in the production run, and the subkeyword USERAWDAT of the key VIBRON in the (evaluation) restart. Such a restart does not invoke any new SCF and will, therefore, typically only take a couple of seconds or minutes. If SAVERAWDAT is not specified, restarts are still possible, but the energy window cannot be adjusted differently, and no new state selection can be performed.

Resonance Raman Input options

A number of options are available for the VIBRON module, most of which are for special applications. All the options mentioned below have to appear in the VIBRON block. The VIBRON module always requires an EXCITATION input block, in which the total number of excited states to be calculated must be specified.

VIBRON
 NMTAPE filename
 RESRAMAN
 {DISPTYPE disptype}
 {STPSIZE stpsize}
 {ONLYSYM}
 {NOTONLYSYM}
 {DOMODES list}
 {DONTMODES list}
 {DSCHEME dscheme}
 {EUTHRES euthres}
 {ELTHRES elthres}
 {SELSTATE list}
 {SAVERAWDAT}
 {USERAWDAT}
END

The most important ones in connection with RR calculations are:

NMTAPE filename
NMTAPE is the only obligatory keyword for the VIBRON module. It specifies the name of a TAPE21 file from a previous frequency calculation. This TAPE21 file is needed to read the normal modes w.r.t. which the derivatives are computed. I.e., a separate frequency calculation must be carried out first.
RESRAMAN
The second subkeyword RESRAMAN invokes the resonance Raman calculation.
DISPTYPE disptype

Select type of displacement steps; possible values are:

  • MASSWE: steps in terms of mass-weighted normal mode vectors [default]
  • CARTES: steps in terms of Cartesian normal mode vectors
  • REDUCE: steps in terms of reduced normal modes
  • ENERGY: steps in terms of expected energy change (according to harmonic approximation)
  • ZPELEV: like energy, but energy is expressed in ZPE units
STPSIZE stpsize
Sets the step size for the numerical differentiation in the default unit for the given DISPTYPE
ONLYSYM
Calculate derivatives only for totally symmetric modes (this is useful since this RR estimate only holds for Franck-Condon type Raman scattering, which is zero for non-symmetric modes). This option is ON per default in RR calculations, it can be switched off with the key NOTONLYSYM.
DOMODES list
Calculate derivatives only for the normal modes with numbers mentioned in list.
DONTMODES list
Calculate derivatives for all normal modes except the ones with numbers mentioned in list.
DSCHEME dscheme

The type of differentiation to be used can be set. Three different values for dscheme are available:

ELCHAR
A simple diabatic picture in which adiabatic states are mapped to the adiabatic states for the references structure based on a maximum transition density overlap criterion [default].
EIGVEC
A diabatic picture in which a diabatization is carried out as explained in Ref. [3]. I.e. the energies used here are the diagonal elements of the potential energy matrix for the nuclear Schrödinger Equation.
ADIABS
The adiabatic picture; can only be used if the symmetry of the excited states is supplied from a separate calculation, since the VIBRON module cannot check which states are allowed to cross as no symmetry is used in the excitation calculations. Consult the ADF-VIBRON manual [4] for information.
ELTHRESH elthresh EUTHRESH euthresh
For advanced users it is possible to set an energy window within the range of states calculated, and only the states within this energy window will be taken into account in the evaluation. - elthresh: lower bound in eV, default 0 eV. - euthresh: upper bound in eV, default 1.0E10 eV.
SELSTATE list
For advanced users it is furthermore possible to pick out certain states from the energy window, and only perform the mapping (and diabatization, if requested) and differentiation for them. Here list includes the number (in ascending excitation energy) of the excited state at the reference (equilibrium) structure.
SAVERAWDAT
All raw data are saved to TAPE21, which might be an enormous amount of data. The selected states and the energy window settings can then be adjusted in a restart (with the usual restart key) and the inclusion of the subkey USERAWDATA after all single-point calculations are done.
USERAWDAT
All raw data are read from a previous calculation. The selected states and the energy window settings can now be adjusted. You need to invoke the usual restart key. Such a restart does not invoke any new SCF and will, therefore, typically only take a couple of seconds or minutes. If SAVERAWDAT is not specified in the previous calculation, restarts are still possible, but the energy window cannot be adjusted differently, and no new state selection can be performed.

VROA: (Resonance) vibrational Raman optical activity

In ADF2010 a method is implemented to calculate both on- and off-resonance vibrational Raman optical activities (VROAs) of molecules using time-dependent density functional theory, see Ref. [5]. This is an extension of a method to calculate the normal VROA by including a finite lifetime of the electronic excited states in all calculated properties. The method is based on a short-time approximation to Raman scattering and is, in the off-resonance case, identical to the standard theory of Placzek. The normal and resonance VROA spectra are calculated from geometric derivatives of the different generalized polarizabilities obtained using linear response theory which includes a damping term to account for the finite lifetime. Gauge-origin independent results for normal VROA have been ensured using either the modified-velocity gauge or gauge-included atomic orbitals. In ADF2016 the velocity gauge tensors required for the calculation of VROA are now correctly calculated with the life time damping parameter. With these complex tensors fixed, resonance VROA intensities are now origin invariant in the velocity gauge, see also Ref. [6].

For the normal VROA use numerical frequencies, and the subkey VROA of the key AORESPONSE. Example input:

GEOMETRY
 Frequencies
 End
END
AORESPONSE
 VROA
 scf  converge 1d-6 iterations 100
 frequency 1 5145 Angstrom
 ALDA
 EL_DIPOLE_EL_DIPOLE VELOCITY
 EL_DIPOLE_EL_QUADRUPOLE VELOCITY
 EL_DIPOLE_MAG_DIPOLE VELOCITY
END

For the resonance VROA use numerical frequencies, and the subkey VROA and LIFETIME of the key AORESPONSE. Example input:

GEOMETRY
 Frequencies
 End
END
AORESPONSE
 VROA
 scf  converge 1d-6 iterations 100
 frequency 1 5.15462 eV
 lifetime 0.0037
 ALDA
 EL_DIPOLE_EL_DIPOLE VELOCITY
 EL_DIPOLE_EL_QUADRUPOLE VELOCITY
 EL_DIPOLE_MAG_DIPOLE VELOCITY
END

In combination with the DIM/QM model one might simulate surface-enhanced Raman optical activity (SEROA), see section on DIM/QM and Ref. [6].

If you use the (deprecated) STO fit option, it’s recommended to include the keyword FitAOderiv in the AORESPONSE block.

Vibrational Circular Dichroism (VCD) spectra

The following keyword enables calculation of rotational strength during an analytical frequencies calculation:

VCD

It is important to note that the VCD keyword only works in combination AnalyticalFreq and symmetry NOSYM. The frozen core approximation and open shell systems are not supported.

AnalyticalFreq
End
SYMMETRY NOSYM
Basis
  type TZP
  core NONE
End

The VCD intensities are calculated using Stephens’ equations for VCD. For the calculation of the atomic axial tensors (AATs), analytical derivatives techniques and London atomic orbitals (the so called GIAO) are employed. As a result the calculated rotational strengths are origin independent, and therefore the common origin gauge is used [7].

Calculation of the AATs requires an analytical frequencies calculation. This limits the choice of functionals that can be used for VCD calculations. See the ANALYTICALFREQ keyword for a complete list of the available functionals. The VCD calculations can be done immediately after a geometry optimization, just like analytical frequencies calculations.

The accuracy of the vibrational rotational strengths are determined by the accuracy of the harmonic force field, atomic polar tensors (APTs) and AATs. The most critical parameter being the harmonic force field. Thus, for a fair comparison with experimental data, accurate geometries and functionals that yield accurate force fields (e.g. BP86, OLYP, etc) should be used. Our tests showed that the BP86 functional in combination with TZP basis sets is always a safe choice. For a comparison of VCD spectra calculated with various functionals (e.g BP86, OLYP, BLYP, B3PW91 and B3LYP) see [7]. Regarding the geometries, we recommend the following strict settings, 10-4 for the geometry convergence of the gradients, and BeckeGrid quality good. The default settings should be used for the calculation of the frequencies.

By default, only the vibrational rotational strengths are printed in the ADF output file. The AATs can also be printed by specifying the keyword:

PRINT VCD

VCDtools is a program that can be used to do an analysis of the VCD spectrum. VCDtools can be used with the ADF-GUI module ADFspectra.

Vibrationally resolved electronic spectra

See the section on vibrationally resolved electronic spectra.

References

[1]L. Jensen, L. Zhao, J. Autschbach and G.C. Schatz, Theory and method for calculating resonance Raman scattering from resonance polarizability derivatives, Journal of Chemical Physics 123, 174110 (2005)
[2]J. Neugebauer, E.J. Baerends, E. Efremov, F. Ariese and C. Gooijer, Combined Theoretical and Experimental Deep-UV Resonance Raman Studies of Substituted Pyrenes, Journal of Physical Chemistry A 109, 2100 (2005)
[3]J. Neugebauer, E.J. Baerends and M. Nooijen, Vibronic coupling and double excitations in linear response time-dependent density functional calculations: Dipole-allowed states of N2 , Journal of Chemical Physics 121, 6155 (2004)
[4]J. Neugebauer, Vibronic Coupling Calculations using ADF, documentation on the VIBRON module available on request.
[5]L. Jensen, J. Autschbach, M. Krykunov, and G.C. Schatz, Resonance vibrational Raman optical activity: A time-dependent density functional theory approach, Journal of Chemical Physics 127, 134101 (2007)
[6](1, 2) D.V. Chulhai and L. Jensen, Simulating Surface-Enhanced Raman Optical Activity Using Atomistic Electrodynamics-Quantum Mechanical Models, Journal of Physical Chemistry A 118, 9069 (2014)
[7](1, 2) V.P. Nicu J. Neugebauer S.K. Wolff and E.J. Baerends, A vibrational circular dichroism implementation within a Slater-type-orbital based density functional framework and its application to hexa- and hepta-helicenes, Theoretical Chemical Accounts 119, 245 (2008)