# Simulating Vibronic Effects in Optical Spectra This tutorial demonstrates how to compute the vibronic spectra of SO$_2$ using CP2K and the VibronicSpec post-processing tool. We'll calculate vibrational modes and excited state forces with CP2K, then combine these results to obtain the vibronic spectrum using different theoretical methods, as shown in the figure. ![](so2_spectrum.png){align=center} The methods developed are based on these publications: - [J. Chem. Phys. 127, 164319 (2007)](https://pubs.aip.org/aip/jcp/article/127/16/164319/937415/Analysis-and-prediction-of-absorption-band-shapes) - [J. Chem. Phys. 137, 234107 (2012)](https://pubs.aip.org/aip/jcp/article/137/23/234107/191739/Efficient-and-automatic-calculation-of-optical) - [J. Chem. Theory Comput. 13, 2823 (2017)](https://pubs.acs.org/doi/10.1021/acs.jctc.7b00325) ## Theoretical Background **Vibronic Transitions and the Franck-Condon Principle** Vibronic spectra arise from transitions between electronic states that involve simultaneous changes in both electronic and vibrational states. The intensity distribution follows the Franck-Condon principle, which states that electronic transitions occur much faster than nuclear motions, leading to vertical transitions where nuclei maintain their initial positions and momenta. **Displaced Harmonic Oscillator Model** Within the harmonic approximation, the potential energy surfaces of ground and excited states are described as parabolas with identical curvatures (same frequencies and normal modes) but displaced minima. The key parameter is the displacement vector $\Delta Q_k$, which represents the geometry change between equilibrium structures along each normal mode coordinate: $$ \Delta Q_k = \frac{1}{\omega_k^2} \frac{\partial E}{\partial Q_k} $$ where $\omega_k$ is the vibrational frequency of mode $k$ and $\partial E/\partial Q_k$ is the energy gradient along the normal mode coordinate. **Huang-Rhys Factors and Relaxation Energy** The Huang-Rhys factor $S_k$ quantifies the vibronic coupling strength for each vibrational mode $k$: $$ S_k = \frac{1}{2}\omega_k\Delta Q_k^2 $$ Large $S_k$ values indicate strong coupling, leading to pronounced vibrational progression in the spectrum. The total relaxation energy $\lambda = \sum_k S_k\omega_k$ represents the energy stabilization due to geometry relaxation after electronic excitation. **Thermal Occupation Factors** The thermal occupation number $\bar{n}_k$ represents the average number of vibrational quanta excited at temperature T: $$ \bar{n}_k = \frac{1}{e^{\omega_k/k_BT} - 1} $$ This factor accounts for thermally excited vibrational states that can participate in the transition, with $\bar{n}_k = 0$ at zero temperature. ### Lineshape Methods The simplest method is to calculate the vertical excitations and convolute each transition with an empirical broadening function. This approach, however, apart from the ambiguity in the choice of the width parameter, does not recognize the asymmetry of the vibronic band. ### IMDHO Method (Independent Mode Displaced Harmonic Oscillator) The IMDHO method uses the time-domain correlation function formalism: $$ \sigma(\omega) = \frac{4\pi^2\omega}{3c} f \int_{-\infty}^{\infty} dt \, e^{i(\omega - E_{ad} - \lambda)t} \prod_{k=1}^{N} \exp\left[-S_k(2\bar{n}_k+1) + S_k\bar{n}_ke^{i\omega_kt} + S_k(\bar{n}_k+1)e^{-i\omega_kt}\right] $$ where: - $f$ is the electronic oscillator strength, - $E_{ad} = E_{vert} - \lambda$ is the adiabatic excitation energy. - The product term describes the vibrational overlap for all modes, - The exponential terms represent creation and annihilation of vibrational quanta, - Temperature effects are included through the $\bar{n}_k$ factors. This method provides the highest accuracy by treating each vibrational mode independently, allowing it to capture full vibronic progression, mode-specific displacements, and temperature effects. ### LQ2 Method - Symmetric Gaussian Model (Second-Order Cumulant Expansion) The LQ2 method is a Gaussian approximation that considers only the second spectral moment: $$ \sigma(\omega) \propto f \exp\left[-\frac{(\omega - E_{vert})^2}{4\alpha}\right] $$ where the linewidth parameter $\alpha$ is defined as: $$ \alpha = \sum_k S_k\omega_k^2 $$ This method provides a symmetric lineshape and is computationally efficient but neglects spectral asymmetry and detailed vibronic structure. ### LQ3 Method - Asymmetric Gaussian Model (Third-Order Cumulant Expansion) The LQ3 method includes anharmonic corrections through a cubic time term: $$ \sigma(\omega) \propto f \int dt \, e^{i(\omega - E_{vert})t} e^{-\alpha t^2 + \gamma t^3} $$ where, in addition to the linewidth parameter $\alpha$ from above, the anharmonicity parameter $\gamma$ then captures spectral asymmetry: $$ \gamma = \frac{1}{12}\sum_k S_k\omega_k^3 $$ This method provides improved accuracy over LQ2 by accounting for spectral skewness while remaining more computationally efficient than the full IMDHO calculation. ### Method Selection The choice of method depends on the desired balance between computational cost and spectral accuracy: - LQ2: Very fast, suitable for preliminary calculations and systems with weak vibronic coupling. - LQ3: Fast, captures asymmetric lineshapes. - IMDHO: Most accurate, essential for detailed vibronic progression and temperature-dependent studies. Can be expensive for large systems and many excited states. All methods utilize the same fundamental physical parameters (frequencies, displacements, Huang-Rhys factors) but employ different levels of approximation in constructing the final spectral lineshape. ## Example: SO$_2$ molecule In this example, we'll use SO$_2$ and calculate the vibronic spectra. To compute the vibronic spectra, we need 2 separate CP2K calculations: - Vibrational analysis of the ground-state - TDDFT absorption spectrum with excited-state forces This tutorial will not provide detailed information regarding these mentioned calculations. Since vibronic spectroscopy methods are highly sensitive to molecular geometry, you must start from a well-converged global minimum. The geometry file provided in this tutorial has been converged with the [MAX_FORCE](#CP2K_INPUT.MOTION.GEO_OPT.MAX_FORCE) keyword set to `1.0E-5` to ensure it meets this requirement. You can check the geometry by looking at the vibrational analysis. If there are negative frequencies, the geometry might not be converged properly. For solids and liquids, however, vibrational analysis may produce negative frequencies that require careful interpretation. Additionally, the electronic structure calculations must use the same parameters, i.e., functional, basis set, and convergence settings, across all calculations. The optimized geometry is specific to the chosen computational setup, so using different parameters in vibrational analysis and spectrum calculation could mean the geometry is no longer at the true minimum for those new settings. The geometry in this tutorial is optimized using `PBE0` functional and `TZVP-MOLOPT-PBE0-GTH` basis sets. The [ADMM](#CP2K_INPUT.FORCE_EVAL.DFT.AUXILIARY_DENSITY_MATRIX_METHOD) approximation ([](#Guidon2010)) is also used with `admm-dzp` basis to reduce the cost of exchange integrals. ### 1. Vibrational modes and frequencies First, we calculate the vibrational spectrum using the following input. SO$_2$ molecule has 3 modes, ($3*N_\textrm{atom}-6$), therefore the [NPROC_REP](#CP2K_INPUT.VIBRATIONAL_ANALYSIS.NPROC_REP) value should be adjusted to the requested processors. Intensities are not required. ```none &GLOBAL PROJECT so2 PRINT_LEVEL LOW RUN_TYPE VIBRATIONAL_ANALYSIS &END GLOBAL &VIBRATIONAL_ANALYSIS INTENSITIES F FULLY_PERIODIC FALSE NPROC_REP 6 &END VIBRATIONAL_ANALYSIS &FORCE_EVAL METHOD Quickstep &DFT CHARGE 0 MULTIPLICITY 1 BASIS_SET_FILE_NAME BASIS_MOLOPT_UZH POTENTIAL_FILE_NAME POTENTIAL_UZH BASIS_SET_FILE_NAME BASIS_ADMM_UZH &QS METHOD GPW EPS_DEFAULT 1.0E-12 MIN_PAIR_LIST_RADIUS -1.0 &END QS &MGRID CUTOFF 600 REL_CUTOFF 60 &END MGRID &SCF SCF_GUESS RESTART EPS_SCF 1.0E-7 MAX_SCF 20 &OT MINIMIZER DIIS PRECONDITIONER FULL_ALL &END &OUTER_SCF EPS_SCF 1.0E-7 MAX_SCF 100 &END &END SCF &XC &XC_FUNCTIONAL PBE0 &END XC_FUNCTIONAL &HF FRACTION 0.25 &HF_INFO ON &END &INTERACTION_POTENTIAL POTENTIAL_TYPE TRUNCATED CUTOFF_RADIUS 3.0 &END &SCREENING EPS_SCHWARZ 1.0E-10 &END &END HF &END XC &AUXILIARY_DENSITY_MATRIX_METHOD METHOD BASIS_PROJECTION ADMM_PURIFICATION_METHOD NONE EXCH_SCALING_MODEL NONE EXCH_CORRECTION_FUNC PBEX &END &POISSON PERIODIC NONE POISSON_SOLVER WAVELET &END POISSON &END DFT &SUBSYS &COORD S 4.9999986448 4.9666767697 5.2571227806 O 4.9999999980 6.2087412933 4.5332461336 O 4.9999999958 3.7245960698 4.5332576512 &END COORD &TOPOLOGY &CENTER_COORDINATES &END &END TOPOLOGY &CELL ABC 10 10 10 ALPHA_BETA_GAMMA 90.0 90.0 90.0 PERIODIC NONE &END CELL &KIND S BASIS_SET TZVP-MOLOPT-PBE0-GTH-q6 BASIS_SET AUX_FIT admm-dzp-q6 POTENTIAL GTH-PBE0-q6 &END KIND &KIND O BASIS_SET TZVP-MOLOPT-PBE0-GTH-q6 BASIS_SET AUX_FIT admm-dzp-q6 POTENTIAL GTH-PBE0-q6 &END KIND &END SUBSYS &END FORCE_EVAL ``` From this calculation, we obtain the vibrational spectrum in `so2-VIBRATIONS-1.mol` in molden format, which we will use for the post-processing. This file contains the atomic coordinates in a.u., frequencies (`[FREQ]`) in $\textrm{cm}^{-1}$ and normal modes (`[FR-NORM-COORD]`) in a.u.. ```none [Molden Format] [Atoms] AU S 1 16 9.448629 9.448646 10.132595 O 2 8 9.448632 11.795808 8.764666 O 3 8 9.448632 7.101454 8.764688 [FREQ] 463.497878 1170.275972 1387.344563 [FR-COORD] S 9.448629 9.448646 10.132595 O 9.448632 11.795808 8.764666 O 9.448632 7.101454 8.764688 [FR-NORM-COORD] vibration 1 -0.000001 0.000001 0.408496 0.000001 -0.499005 -0.409341 0.000001 0.499004 -0.409341 vibration 2 -0.000001 0.000037 0.344976 0.000001 0.566526 -0.345673 0.000001 -0.566601 -0.345707 vibration 3 -0.000000 -0.520552 0.000027 0.000001 0.521672 -0.304029 -0.000001 0.521585 0.303974 ``` ### 2. Absorption spectrum and excited-state forces The absorption spectrum and excited‑state forces can be calculated in a single run. For this, we will change the [RUN_TYPE](#CP2K_INPUT.GLOBAL.RUN_TYPE) to `ENERGY_FORCE`, remove the `VIBRATIONAL_ANALYSIS` section from the input above, and add the [PROPERTIES](#CP2K_INPUT.FORCE_EVAL.PROPERTIES) section below to `FORCE_EVAL`. ```none &PROPERTIES &TDDFPT NSTATES 10 CONVERGENCE [eV] 1.0e-5 ADMM_KERNEL_CORRECTION_SYMMETRIC &PRINT &FORCES THRESHOLD 0.01 &END &END PRINT &END TDDFPT &END PROPERTIES ``` For this example, we’ll request 10 excited states, which captures the energy range relevant to our spectrum. For other molecules, this number should be adjusted based on the system size and the energy range of interest. We also set the excitation convergence to $10^{-5}$ eV, making it tighter than the default, to ensure more accurate excitation energies. We also ask to print the forces, and we set a [THRESHOLD](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.TDDFPT.PRINT.FORCES.THRESHOLD) of 0.001. This means forces are only calculated for states whose oscillator strength is above that value. Filtering out weak states helps keep the calculation efficient. You can also list specific states explicitly with the [LIST](#CP2K_INPUT.FORCE_EVAL.PROPERTIES.TDDFPT.PRINT.FORCES.LIST) keyword ([Strand2019](https://doi.org/10.1063/1.5078682), [](#Iannuzzi2005), [](#Hehn2022)). The output file we need is `so2-TDFORCE-1_0.tdfrc` file, which contains the spectrum data and the forces in this format: ```none # STATE NR. 1 4.42928 eV Oscillator strength: 0.0018782 3 -0.00000011533194 -0.00001147372289 0.08283325195246 -0.00000003787961 0.08521318053390 -0.03995448825912 -0.00000002363537 -0.08518437233452 -0.04121974562195 # STATE NR. 2 4.73217 eV Oscillator strength: 0.0000000 # STATE NR. 3 6.88229 eV Oscillator strength: 0.0266649 3 -0.00000029284682 -0.00000325607886 0.16080039585134 0.00000007170914 0.09720148772146 -0.07996409950757 0.00000007474609 -0.09720032431716 -0.07989742343589 # STATE NR. 4 8.59178 eV Oscillator strength: 0.0000123 # STATE NR. 5 8.82238 eV Oscillator strength: 0.0000000 # STATE NR. 6 8.96389 eV Oscillator strength: 0.0508470 3 -0.00000000259255 -0.00000738331081 -0.01039580009444 -0.00000007113333 0.13252622525927 0.00524117205768 -0.00000007174747 -0.13251465091338 0.00488212051205 # STATE NR. 7 9.69782 eV Oscillator strength: 0.0000220 # STATE NR. 8 9.73440 eV Oscillator strength: 0.0066963 3 -0.00000037994327 -0.00002119798164 0.20468077319361 0.00000011356579 0.10583044117747 -0.10298557213750 0.00000011155504 -0.10580720493312 -0.10318459031726 # STATE NR. 9 9.77699 eV Oscillator strength: 0.0000000 # STATE NR. 10 10.99706 eV Oscillator strength: 0.0150407 3 -0.00000061611618 0.00000032773220 0.30829793267287 0.00000029373336 0.08008916703722 -0.15453299977376 0.00000029256088 -0.08010428717480 -0.15324331957030 ``` The forces for the states with an intensity higher than the cutoff are calculated and printed, whereas less intense states are skipped. ### 3. Calculating vibronic spectrum We will use the python package [VibronicSpec](https://github.com/cp2k/cp2k/tree/master/tools) that can be found in CP2K tools. An example config.toml file is provided. We will modify it for the SO$_2$ molecule. ```none [files] vibrations_filename = "so2-VIBRATIONS-1.mol" output_filename = "so2_spectrum.txt" force_filename = "so2-TDFORCE-1_0.tdfrc" [calculation] # Single key. Choose one format: # states = "all" # Use all available states with forces # states = [1, 2, 6] # Explicit state list states = "threshold:0.001" # Filter oscillator strength by threshold energy_min = 4 # in eV energy_max = 11 # in eV energy_points = 1000 method = "lq2" spectrum_type = "absorption" stokes_shift = 0.0 [output] print_individual_states = false [imdho] temperature = 298 # in K gamma_broadening = 0.012398418754199978 # in eV theta_broadening = 0.0 [integration] max_time_slices = 5000 lq3_time_step = 30 lq3_convergence = 0.0000001 imdho_time_step = 30 imdho_convergence = 0.000000001 cutoff_energy = 0.0816 ``` The CP2K output files from the vibrational analysis and the TDDFT forces should be provided. Key calculation parameters are: - Energy range (in eV) - Number of points for the spectrum - Method (LQ2, LQ3, or IMDHO) - Spectrum type (absorption or fluorescence) - States to include, it can take different arguments: - "all", includes all states that have forces in TDFORCE file - a list like [1, 2, 5] to explicitly give state numbers - "threshold:0.001", this option works like the THRESHOLD filter we used in the TDDFT force calculation, including only states with oscillator strength above 0.001. Optional sections provide defaults for method‑specific settings. For the IMDHO method, you can set the temperature (in K) for thermal occupation factors and also the gamma broadening (in eV). Integration parameters control the time integrations required for LQ3 and IMDHO. The defaults are tuned values and they should be used unless you know what you’re doing. The final output is a spectrum summed over all included states. If you want to see each state separately, set `print_individual_states = true`. VibronicSpec will print information regarding the calculation, such as the states filtered, which should be checked. It will also print the progression, and total runtime. The figure above shows the calculated absorption spectrum of SO$_2$ using all three methods using the provided input, only changing the method. The LQ2 methods gives a simple Gaussian profile without the vibronic structure. The LQ3 methods adds a cubic coupling through time integration, revealing some shoulder features. The IMDHO method (calculated at 298 K) includes full mode-by-mode displacement and temperature effects, producing the most detailed spectrum. The LQ2 and LQ3 intensities are scaled to align with the IMDHO peak for visual comparison. For reference, for this particular example, if LQ2 takes one unit of time, then LQ3 and IMDHO will take about 800× longer for this molecule with the given setup. In practice, timings scale with the number of excited states, vibrational modes, and spectral resolution.