UV-Vis Spectra from ORCA

5 minute read

Published:

I’ve been using TD-DFT to simulate some UV-Vis spectra for a collaborator, and wanted to understand what ORCA’s orca_mapspc utility is actually doing.Here’s how to plot exactly the same spectrum as orca_mapspc.

Running the calculation

The ORCA manual, and the ORCA input library should be consulted for instructions on how to run a TD-DFT calculation, or any sort of CI calculation.

Understanding the result

ORCA will print the information you need to the output file in a section with the following header

         ABSORPTION SPECTRUM VIA TRANSITION ELECTRIC DIPOLE MOMENTS
-----------------------------------------------------------------------------
State   Energy    Wavelength  fosc         T2        TX        TY        TZ
        (cm-1)      (nm)                 (au**2)    (au)      (au)      (au)
-----------------------------------------------------------------------------

We only need the Energy and fosc entries to plot a spectrum. The fosc column contains the oscillator strength of each transition. In order to understand what this means, recall the integrated absorption coefficient of a given 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}]$.

The oscillator strength $f$ is defined as the ratio of $\mathcal{A}_\mathrm{e}$ for a given band to the (constant) value for that of a harmonically oscillating electron.

\[f = \frac{\int_\mathrm{band}{\epsilon_\mathrm{e}\left(\widetilde{\nu}\right)\ \mathrm{d}\widetilde{\nu}}}{\mathcal{A}_\mathrm{e}^{\mathrm{electron}}} = \frac{\mathcal{A}_\mathrm{e}^{\mathrm{ORCA}}}{\mathcal{A}_\mathrm{e}^{\mathrm{electron}}}\]

and usually has a value between 0 and 1. In order to plot the spectrum, we need to convert $f$ back into $\mathcal{A}_\mathrm{e}^\mathrm{ORCA}$ since this tells us how large each absorption band will be.

\[\mathcal{A}_\mathrm{e}^\mathrm{ORCA} \ [\mathrm{1000 \ cm \ mol^{-1}}] = f \cdot \mathcal{A}_\mathrm{e}^{\mathrm{electron}} = f \cdot 2.31\times 10^8 \ [\mathrm{1000 \ cm \ mol^{-1}}]\]

An explanation of the oscillator strength can be found on Pages 80 and 81 of Barrow’s Introduction to Molecular Spectroscopy.

Plotting a spectrum

To generate the UV-Vis spectrum orca_mapspc uses an area-normalised Lorentzian lineshape $L(x)$ for each transition

\[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 transition.

For a given transition, the corresponing Lorentzian is multiplied by the linear integrated absorption coefficent $\mathcal{A}^\mathrm{ORCA}$ obtained by dividing $\mathcal{A}_\mathrm{e}^\mathrm{ORCA} \ [\mathrm{1000 \ cm \ mol^{-1}}]$ by $\ln(10)$. The multiplication of $L(x) \ [\mathrm{cm}]$ by $\mathcal{A}^\mathrm{ORCA} \ \ [\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 transitions as a single spectrum, simply sum $L(x) \cdot \mathcal{A}^\mathrm{ORCA}$ for each transition, and plot against wavenumber.

What about Velocity?

You might have noticed that ORCA prints an identical absorption spectrum table with the title

         ABSORPTION SPECTRUM VIA TRANSITION VELOCITY DIPOLE MOMENTS

which contains an identical set of transitions with slightly different oscillator strengths. This stems from the use of a different operator in calculating transition intensities which should give equivalent results, but rarely does due to the approximations used in quantum chemistry.

In order to calculate the oscillator strength of a transition between two states we need to know the value of the transition dipole moment for that transition.

The transition dipole moment is a matrix element of the dipole operator

\[\mu_{ba} = \langle \psi_b | \hat{\mu} | \psi_a \rangle\]

In the case of the electric dipole moment and a particle in three-dimensions

\[\mu_{ba} = \langle \psi_b | q\mathbf{r} | \psi_a \rangle = q\int \psi_b^*(r) \ \mathbf{r} \ \psi_a(\mathbf{r}) \ \mathrm{d}r\]

where $q$ is the charge of the particle, and $\mathbf{r}$ the position operator.

The position operator and Hamiltonian operator have the well-known commutation relation

\[\left[\hat{H}, \mathbf{r}\right] = -\frac{\hbar}{m}\vec{\nabla}\]

Which can be obtained from the derivation for a single dimension \(\left[\hat{H}, x\right] = \left[\frac{p_x^2}{2m} + V, x\right] = \left[\frac{p_x^2}{2m}, x\right] + \left[V, x\right] = \frac{\hbar}{im} p_x\)

where \(\left[\frac{p_x^2}{2m}, x\right] = \frac{1}{2m}\left\{p_x[p_x,x] + [p_x,x]p_x\right\}\) \([p_x,x] = -i\hbar\)

Alternatively $\left[\hat{H}, \mathbf{r}\right]$ can be evaluated directly

\[\left \langle a \left | \left[\hat{H}, \mathbf{r}\right] \right | b \right \rangle = \left \langle a \left | \hat{H} \mathbf{r} \right | b \right \rangle - \left \langle a \left | \mathbf{r} \hat{H} \right | b \right \rangle = \left(E_a = E_b \right) \left \langle a \left | \mathbf{r} \right | b \right \rangle\]

so then

\[\left(E_a = E_b \right) \left \langle a | \mathbf{r} | b \right \rangle = -\frac{\hbar}{m}\left \langle a \left | \vec{\nabla} \right | b \right \rangle\]

Meaning that one could obtain transition dipole moment values by either computing the position operator directly, or by computing the velocity operator and using the above relationship.

Due to the number of approximations used in computational chemistry these two approaches rarely give the same answer, and so the user is free to choose which one they want to use. orca_mapspc allows the user to specifiy which form to use, but by default uses the position values.