UV-Vis Spectra from ORCA
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
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
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\}\) and \([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 we can 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. Therefore, ORCA and orca_mapspc allow the user to choose which value to use, and by default orca_mapspc uses the position form which they refer to as the electric dipole.