# Resonance Raman (RR) and absorption spectra `vibrant` can compute static or molecular dynamics (MD)-based resonance Raman (RR) and absorption spectra processing the induced dipole moments calculated via real-time time-dependent density functional theory (RT-TDDFT). [RT-TDDFT](https://pubs.acs.org/doi/10.1021/acs.chemrev.0c00223) models resonance effects by propagating the electron density in real time under a small external perturbation, enabling the full excitation spectrum to be captured from a single simulation via Fourier transformation of the time-dependent polarizabilities. This is an important advantage when compared to the more commonly used linear-response TDDFT (LR-TDDFT), which requires specifying a subset of excitations. As a result, the RT-TDDFT-based approach is efficient and easily combinable with both static and MD-based simulations. `vibrant` calculates the RR intensities for a range of incident laser wavelengths $\tilde{\nu}_{in}$ specified in the input file via the `laser_in` keyword in the `raman` section. Every RR spectrum calculated for a different laser wavenumber are appended as additional columns. More details on the differences of static and MD-based RR and absorption spectra are provided in the following sections, together with the details on the Padé post-processing which significantly enhances the spectral resolution. ## a) Static RR and absorption intensities The static RR calculation in `vibrant` relies on the normal mode analysis within the harmonic approximation. The static RR intensities are computed from the derivatives of the frequency-dependent polarizability tensors $\boldsymbol{\alpha}$ with respect to the mass-weighted normal coordinates $Q_{p}$. Following [Placzek's polarization theory](https://onlinelibrary.wiley.com/doi/book/10.1002/0470845767), the [RR intensity](https://pubs.acs.org/doi/full/10.1021/acs.jctc.9b00512) at laser frequency $\omega$ for each normal mode frequency $\tilde{\nu}_{p}$ is given as: $$ I_{\textnormal{RR}} (\tilde{\nu}_{p}, \omega) = \frac{h}{8\epsilon_{0}^{2}c} \frac{(\omega-\tilde{\nu}_{p})^{4}}{\tilde{\nu}_{p} \left(1-\exp\left(-\frac{hc\tilde{\nu}_{p}}{k_{B}T}\right)\right)} \frac{X\beta^{2}_{p} (\omega)+Y\gamma_{p}^{2}(\omega) }{45} $$ where $h$ is the Planck’s constant, $\varepsilon_{0}$ is the vacuum permittivity, $c$ is speed of the light, $k_{B}$ is the Boltzmann constant, $T$ is the temperature. [$X$ and $Y$](https://onlinelibrary.wiley.com/doi/book/10.1002/0470845767) can take different values depending on the polarization configuration of the scattered light as explained in Section [Raman spectra](Raman_spec.md). $\beta^{2}_{p}(\omega)$ and $\gamma_{p}^{2}(\omega)$ correspond to the frequency-dependent isotropic and anisotropic contributions, respectively. They are computed from the derivatives of the dynamic polarizability tensor $\boldsymbol{\alpha}(\omega ,\tilde{\nu}_p)$ with respect to the normal mode coordinate $Q_{p}$. $\beta^{2}_{p}(\omega)$ includes the diagonal components of $\boldsymbol{\alpha}(\omega ,\tilde{\nu}_p)$: $$ \beta^{2}_{p}(\omega)=\frac{1}{9}\left ( \frac{\partial \alpha _{xx}(\omega)}{\partial Q_{p}}+\frac{\partial\alpha _{yy}(\omega) }{\partial Q_{p}} +\frac{\partial \alpha _{zz}(\omega)}{\partial Q_{p}}\right )^{2} $$ And $\gamma_{p}^{2}(\omega)$ includes also the off-diagonal terms: $$ \gamma_{p}^{2}(\omega) = \frac{1}{2}\left(\frac{\partial \alpha_{xx}}{\partial Q_{p}} - \frac{\partial \alpha_{yy}}{\partial Q_{p}}\right)^{2} + \frac{1}{2}\left(\frac{\partial \alpha_{yy}}{\partial Q_{p}} - \frac{\partial \alpha_{zz}}{\partial Q_{p}}\right)^{2} + \frac{1}{2}\left(\frac{\partial \alpha_{zz}}{\partial Q_{p}} - \frac{\partial \alpha_{xx}}{\partial Q_{p}}\right)^{2} + 3\left(\frac{\partial \alpha_{xy}}{\partial Q_{p}}\right)^{2} + 3\left(\frac{\partial \alpha_{yz}}{\partial Q_{p}}\right)^{2} + 3\left(\frac{\partial \alpha_{zx}}{\partial Q_{p}}\right)^{2} $$ The [absorption intensities](https://pubs.aip.org/aip/jcp/article/149/17/174108/197213) are computed from the trace of the frequency-dependent absorption cross-section tensor $\boldsymbol{\sigma}(\omega)$: $$ I_{\textnormal{abs}} (\omega) = \frac{1}{3}\textnormal{Tr}\left[ \boldsymbol{\sigma}(\omega) \right] $$ where $\boldsymbol{\sigma}(\omega)$ contains only the diagonal components of the imaginary part of the dynamic polarizability tensor $\alpha(\omega)$, as given by: $$ \sigma_{ii} (\omega) = \frac{4 \pi \omega}{c}\textnormal{Im}\left[ \alpha_{ii} (\omega) \right] \quad \textnormal{with} \quad i = x, y, z $$ The dynamic polarizability tensor $\boldsymbol{\alpha}(\omega)$ is computed from the frequency-dependent induced dipole moment ${\mu}^{j}_{i}(\omega)$: $$ \alpha_{ij}(\omega) = \frac{{\mu}^{j}_{i}(\omega) }{E} \quad \textnormal{with} \quad i, j = x, y, z $$ where ${E}$ is an electric field applied along the $x$-, $y$-, or $z$-direction. ${\mu}^{j}_{i}(\omega)$ is obtained via applying a fast Fourier transform to the time-dependent induced dipoles using the [FFTW](https://www.fftw.org/) library: $$ {\mu}^{j}_i(\omega) = \int_{0}^{T} \left( {\mu}^{j}_i(t) - \mu^{0}_{i} \right) e^{i \omega t} e^{- \Gamma t} \mathrm{d}t $$ where $T$ is the total RT-TDDFT simulation time, $t$ is the RT-TDDFT time step, ${\mu}^{j}_i(t)$ is the time-dependent dipole moment under the electric field applied in the $j$ direction, $\mu^{0}_{i}$ denotes the initial static dipole moment of the unperturbed system, and $\Gamma$ is the [damping factor](https://pubs.aip.org/aip/jcp/article/138/4/044101/193408). To compute a static RR spectrum, the user must provide the time-dependent dipole moments for each displaced structure, combined into a single file (see [File Formats](file_formats.md#62-static-spectra-based-on-normal-modes) for details). There should be three files in total, each corresponding to calculations performed under electric fields applied along the x-, y-, and z-directions. An example input section for performing a static RR calculation is shown below: ```bash &global spectra RR temperature fwhm &end global &system filename &end system &static ... &end static &dipoles type_dipole berry field_strength dip_x_file dip_y_file dip_z_file &end dipoles &rtp rtp_time_step rtp_framecount damping_constant &end rtp ... ``` The `rtp` section provides the necessary information about the RT-TDDFT parameters. The `static` section should be given similar to the [static Raman](Raman_spec.md#a-static-raman-intensities) input. Complete input files can be found at [Examples](Examples.md). `vibrant` applies Gaussian broadening to the final discrete set of frequencies and intensities. The `fwhm` keyword controls the full width at half-maximum (FWHM) value in cm ${^{-1}}$ and if not specified, it is set to 10 cm ${^{-1}}$. ```{note} The keyword `write_mol_file` is optional and it executes the printing of a `.mol` file, which includes the optimized geometry (Å), normal mode frequencies (cm$^{-1}$), the non-mass-weighted eigenvectors in atomic units (Bohr) and the non-broadened Raman intensities. The `.mol` file can be opened with [MOLDEN](https://www.theochem.ru.nl/molden/) to visualize the normal modes alongside the Raman spectrum. If multiple incident laser frequencies are requested, `vibrant` generates a separate `.mol` file for each frequency. ``` The absorption spectrum calculation does not require the `static` section. An example input section is given below: ```bash &global spectra ABS &end global &system ... &end system &polarizabilities type_pol induced field_strength &end polarizabilities &dipoles type_dipole berry dip_x_file dip_y_file dip_z_file &end dipoles &rtp rtp_time_step rtp_framecount damping_constant &end rtp ``` The final static resonance Raman intensities are reported in 10 $^{-30}$ cm $^2$ per molecule/system. The absorption spectrum frequencies are reported in eV and the intensities are reported in atomic units. ## b) MD-based RR intensities In an MD-based RR calculation, the polarizability derivatives along the normal mode coordinates $Q_{p}$ are replaced with the [time derivatives of the polarizability autocorrelation functions](https://pubs.rsc.org/en/content/articlehtml/2013/cp/c3cp44302g) similarly to the MD-based Raman calculation. We apply fast Fourier transformations to the time-domain dynamic polarizability autocorrelations using the [FFTW](https://www.fftw.org/) library to convert them into the frequency-domain. The [MD-based RR intensities](https://pubs.acs.org/doi/full/10.1021/acs.jctc.9b00512) at laser frequency $\omega$ is expressed as a linear combination of isotropic and anisotropic polarizabilities: $$ I_{\textnormal{RR}} (\tilde{\nu}, \omega) = \frac{2h}{8\epsilon_{0}^{2}k_{B}T} \frac{(\omega-\tilde{\nu})^{4}}{\tilde{\nu}} \frac{1}{1-\exp\left(-\frac{hc\tilde{\nu}}{k_{B}T}\right)}\frac{X\delta^{2} (\tilde{\nu}, \omega)+Y\epsilon^{2} (\tilde{\nu}, \omega) }{45} $$ where $h$ is the Planck’s constant, $\varepsilon_{0}$ is the vacuum permittivity, $k_{B}$ is the Boltzmann constant, $c$ is the speed of light and $T$ is the temperature which is specified by the user with the keyword `temperature` in `global` section. $X$ and $Y$ are functions of the scattering geometry and polarization as discussed in the section above. Similarly to the static approach, $\delta^{2} (\tilde{\nu})$ and $\epsilon^{2} (\tilde{\nu})$ refer to the [isotropic and anisotropic contributions](https://pubs.acs.org/doi/full/10.1021/acs.jctc.9b00512): $$ \begin{aligned} \delta^2(\tilde{\nu}, \omega) &= \int_{0}^{\hat{T}} \left\langle \frac{\dot{\alpha}_{xx}(\omega, \tau) + \dot{\alpha}_{yy}(\omega, \tau) + \dot{\alpha}_{zz}(\omega, \tau)}{3} \cdot \frac{\dot{\alpha}_{xx}(\omega, \tau + t) + \dot{\alpha}_{yy}(\omega, \tau + t) + \dot{\alpha}_{zz}(\omega, \tau + t)}{3} \right\rangle_{\tau} \, e^{-2\pi i c \tilde{\nu} t} \, dt \end{aligned} $$ where $\tau$ and $t$ refer to the times for two different MD snapshots and $\hat{T}$ is the total MD simulation time. $\dot{\boldsymbol{\alpha}}$ denotes the time derivative of the polarizability tensor $\boldsymbol{\alpha}$. The anisotropic contribution $\epsilon^{2} (\tilde{\nu})$ is given by: $$ \begin{aligned} \epsilon^²(\tilde{\nu}, \omega) = \int_{-\infty}^{\infty} \Biggl[ &\tfrac{1}{2} \left\langle \{\dot{\alpha}_{xx}(\omega, \tau) - \dot{\alpha}_{yy}(\omega, \tau)\} \{\dot{\alpha}_{xx}(\omega, \tau + t) - \dot{\alpha}_{yy}(\omega, \tau + t)\} \right\rangle_{\tau} \\[4pt] &+ \tfrac{1}{2} \left\langle \{\dot{\alpha}_{yy}(\omega, \tau) - \dot{\alpha}_{zz}(\omega, \tau)\} \{\dot{\alpha}_{yy}(\omega, \tau + t) - \dot{\alpha}_{zz}(\omega, \tau + t)\} \right\rangle_{\tau} \\[4pt] &+ \tfrac{1}{2} \left\langle \{\dot{\alpha}_{zz}(\omega, \tau) - \dot{\alpha}_{xx}(\omega, \tau)\} \{\dot{\alpha}_{zz}(\omega, \tau + t) - \dot{\alpha}_{xx}(\omega, \tau + t)\} \right\rangle_{\tau} \\[4pt] &+ \tfrac{3}{4} \left\langle \{\dot{\alpha}_{xy}(\omega, \tau) + \dot{\alpha}_{yx}(\omega, \tau)\} \{\dot{\alpha}_{xy}(\omega, \tau + t) + \dot{\alpha}_{yx}(\omega, \tau + t)\} \right\rangle_{\tau} \\[4pt] &+ \tfrac{3}{4} \left\langle \{\dot{\alpha}_{yz}(\omega, \tau) + \dot{\alpha}_{zy}(\omega, \tau)\} \{\dot{\alpha}_{yz}(\omega, \tau + t) + \dot{\alpha}_{zy}(\omega, \tau + t)\} \right\rangle_{\tau} \\[4pt] &+ \tfrac{3}{4} \left\langle \{\dot{\alpha}_{zx}(\omega, \tau) + \dot{\alpha}_{xz}(\omega, \tau)\} \{\dot{\alpha}_{zx}(\omega, \tau + t) + \dot{\alpha}_{xz}(\omega, \tau + t)\} \right\rangle_{\tau} \Biggr] \\[4pt] &\times e^{-2\pi i c \tilde{\nu} t} \, dt \end{aligned} $$ Since the dynamic polarizability tensor $\boldsymbol{\alpha}(\omega, t)$ is complex-valued (note that $t$ indicates here the MD time step), the above-autocorrelation functions require the handling of [complex autocorrelations](https://pubs.acs.org/doi/full/10.1021/acs.jctc.9b00512) of any function $f(t)$ as shown below: $$ \begin{aligned} \left\langle f(t) \cdot f(t+\tau) \right\rangle_{\tau} &= \left\langle \mathrm{Re}\!\left(f(t)\right) \cdot \mathrm{Re}\!\left(f(t+\tau)\right) \right\rangle_{\tau} + \left\langle \mathrm{Im}\!\left(f(t)\right) \cdot \mathrm{Im}\!\left(f(t+\tau)\right) \right\rangle_{\tau} \\[6pt] &\quad +\, i \Big[ \left\langle \mathrm{Re}\!\left(f(t)\right) \cdot \mathrm{Im}\!\left(f(t+\tau)\right) \right\rangle_{\tau} - \left\langle \mathrm{Im}\!\left(f(t)\right) \cdot \mathrm{Re}\!\left(f(t+\tau)\right) \right\rangle_{\tau} \Big] \end{aligned} $$ $\boldsymbol{\alpha}(\omega, t)$ is computed from the frequency-dependent induced dipole moment ${\mu}^{j}_{i}(\omega, t)$: $$ \alpha_{ij}(\omega, t) = \frac{{\mu}^{j}_{i}(\omega, t) }{E} \quad \textnormal{with} \quad i, j = x, y, z $$ where ${E}$ is an electric field applied along the $x$-, $y$-, or $z$-direction. ${\mu}^{j}_{i}(\omega, t)$ is obtained via applying a fast Fourier transform to the time-dependent induced dipoles using the [FFTW](https://www.fftw.org/) library: $$ {\mu}^{j}_i(\omega, t) = \int_{0}^{T} \left( {\mu}^{j}_i(\tau, t) - \mu^{0}_{i}( t) \right) e^{i \omega \tau} e^{- \Gamma \tau} \mathrm{d}\tau $$ where $T$ is the total RT-TDDFT simulation time, $\tau$ is the RT-TDDFT time step, $t$ is the MD time step, ${\mu}^{j}_i(\tau, t)$ is the time-dependent dipole moment under the electric field applied in the $j$ direction, $\mu^{0}_{i}(t)$ denotes the initial static dipole moment of the unperturbed system, and $\Gamma$ is the damping factor. The MD-based RR calculation in `vibrant` is not fully tested yet, and more tests are being performed. The current implementation requires that the user provides the time-dependent dipole moments for each MD snapshot, combined into a single file (see [File Formats](file_formats.md#61-md-based-spectra) for more details). There should be three files in total, each corresponding to calculations performed under electric fields applied along the x-, y-, and z-directions. An example input section for performing an MD-based RR calculation is shown below: ```bash &global spectra MD-RR temperature &end global &polarizabilities type_pol induced field_strength &end polarizabilities &dipoles type_dipole berry dip_x_file dip_y_file dip_z_file &end dipoles &md time_step correlation_depth &end md &rtp rtp_time_step rtp_framecount damping_constant &end rtp ``` The MD-based absorption spectrum is computed in the same manner as the static absorption spectrum described above, except that the absorption intensities are averaged over all MD snapshots. The MD-based absorption spectrum does not require the specification of the `md` section. ```{note} `vibrant` applies the same post-processing to the final MD-based RR intensities as given in the subsection [Spectral refinement](frequency.md#spectral-refinement), including the application of data mirroring and Hann Window function to the autocorrelation data and the application of the sinc function to the final intensities. ``` The final MD-based RR intensities are reported in m $^2$ K cm 10 $^{-30}$. ```{note} The `write_acf_file` keyword can optionally be used in the `md` section and it executes the printing of a file which contains the time-autocorrelation data. In the case of RR spectrum calculation, this would be the real and imaginary components of the isotropic and anisotropic polarizability-autocorrelation data. ``` ## c) Application of Padé interpolation The spectral resolution of the dynamic polarizability $\alpha(\omega)$ and the resulting absorption spectrum depends on the RT-TDDFT simulation length, which requires hundreds of thousands of time steps for fine energy resolution. To overcome the significant computational cost, `vibrant` can optionally use rational Padé approximants to reconstruct $\alpha(\omega)$ analytically from discrete RT-TDDFT data. The complex-valued $\alpha(\omega)$ is approximated as: $$ \alpha(\omega) \approx \frac{A_0 + A_1 \omega + A_2 \omega^2 + \cdots + A_{\frac{T-1}{2}} \omega^{\frac{T-1}{2}}} {1 + B_1 \omega + B_2 \omega^2 + \cdots + B_{\frac{T}{2}} \omega^{\frac{T}{2}}} $$ where $A_i$ and $B_i$ are the Padé coefficients, and $T$ represents the total number of RT-TDDFT time steps. These coefficients are obtained using Thiele’s reciprocal difference method, applied to the $T$ discrete data points ${\omega_i, \alpha_i}$ generated along the RT-TDDFT trajectory. Padé interpolation improves the resolution and stability significantly for shorter or noisy RT-TDDFT trajectories. Unlike [previous approaches](https://pubs.acs.org/doi/full/10.1021/acs.jctc.6b00511) that replace FFTs entirely, this method applies the Padé interpolation as a post-processing step to refine spectra and correct distortions and noise in the peaks. The Padé fitting is implemented using the "plain-128" algorithm of the [GreenX library](https://nomad-coe.github.io/greenX/gx_ac.html), which provides a numerically robust and efficient solver for [computing Padé coefficients and evaluating the resulting analytic continuation](https://joss.theoj.org/papers/10.21105/joss.07859). The following figure demonstrates how the application of Padé approximants enhances the spectral resolution significantly and "corrects" the peak positions:

Absorption spectra of naphthalene obtained from RT-TDDFT trajectories of varying simulation lengths: (a) without Padé approximants and (b) with Padé approximants.
The upper figure (a) shows the static absorption spectra of naphthalene obtained from RT-TDDFT trajectories of different lengths, where shorter simulations (e.g., 24.2 fs) exhibit peak shifts toward higher energies due to limited spectral resolution. The lower figure (b) demonstrates that applying Padé approximants restores the correct peak positions and yields a converged spectrum even for the shortest trajectory, reducing the total simulation time by roughly a factor of five. The user can request the Padé post-processing for the static or MD-based RR or absorption spectra using the `check_pade` keyword as shown below: ```bash .... &rtp ... check_pade y pade_framecount ... &end rtp ``` The user must also specify the number of the final data points after the Padé fitting using the keyword `pade_framecount`. This parameter controls the resolution of the interpolated data. ```{tip} The Padé fitting usually takes time, but it is accelarated using OpenMP threads. To control the number of parallel threads you can export the environment variable `export OMP_NUM_THREADS=` before calling the `vibrant` code. Occasionally for large RTP steps you can run into a "out-of-memory" error. The amout of RAM memory `vibrant` needs can be computed roughly by $N_{RTP}^2\cdot 16 \cdot 10^{-9} \cdot N_{threads}$ (in GB). ``` More information on the all available keywords can be found on [Keyword Glossary](Keyword_Glossary.rst) and all complete example input files are available on [Examples](Examples.md).