On-the-fly ab initio semiclassical evaluation of absorption spectra of polyatomic molecules beyond the Condon approximation

TL;DR

This paper extends the on-the-fly ab initio thawed Gaussian approximation to include Herzberg-Teller effects, enabling more accurate vibronic spectra calculations beyond the Condon approximation for polyatomic molecules.

Contribution

It introduces a method that incorporates Herzberg-Teller contributions into the thawed Gaussian approximation, improving spectral accuracy beyond standard harmonic models.

Findings

Calculated spectra match experimental data more closely.

Anharmonicity is more significant than Herzberg-Teller in phenyl radical.

Herzberg-Teller effects are crucial for benzene's forbidden transition.

Abstract

To evaluate vibronic spectra beyond the Condon approximation, we extend the on-the-fly ab initio thawed Gaussian approximation by considering the Herzberg-Teller contribution due to the dependence of the electronic transition dipole moment on nuclear coordinates. The extended thawed Gaussian approximation is tested on electronic absorption spectra of phenyl radical and benzene: Calculated spectra reproduce experimental data and are much more accurate than standard global harmonic approaches, confirming the significance of anharmonicity. Moreover, the extended method provides a tool to quantify the Herzberg-Teller contribution: we show that in phenyl radical, anharmonicity outweighs the Herzberg-Teller contribution, whereas in benzene, the Herzberg-Teller contribution is essential, since the transition is electronically forbidden and Condon approximation yields a zero spectrum. Surprisingly, both adiabatic harmonic spectra outperform those of the vertical harmonic model, which describes the Franck-Condon region better. Finally, we provide a simple recipe for orientationally averaging spectra, valid beyond Condon approximation, and a relation among the transition dipole, its gradient, and nonadiabatic coupling vectors.

Equations

$${\sigma }_{\mu \mu }(\omega )=\frac{2\pi \omega }{\hslash c}{\int }_{-\infty }^{\infty }{C}_{\mu \mu }(t)e^{i\omega t}dt$$

$$C_{\mu \mu }(t)=⟨1,g|e^{i{\hat{H}}_{1}t/\hslash }\hat{\mu }{e}^{-i{\hat{H}}_{2}t/\hslash }\hat{\mu }|1,g⟩=⟨\varphi (0)|\varphi (t)⟩e^{-i{E}_{1,g}t/\hslash },$$

$$\overline{\sigma (\omega )}=\frac{1}{3}\text{Tr}[{\sigma }_{\vec{\mu }\vec{\mu }}(\omega )]=\frac{1}{3}\left({\sigma }_{{\mu }_x{\mu }_x}+{\sigma }_{{\mu }_y{\mu }_y}+{\sigma }_{{\mu }_z{\mu }_z}\right).$$

$$\overline{{\sigma }_{\text{FC}}(\omega )}=\frac{1}{3}{\sigma }_{|\vec{\mu }||\vec{\mu }|}(\omega ).$$

$$\vec{\mu }(q)\approx \vec{\mu }(q_0)+{\text{grad}}_q\vec{\mu }{|}_{q_0}^{\top }\cdot (q-q_0)$$

$$\frac{\partial }{\partial q_i}{\vec{\mu }}_{\alpha \beta }=\sum _{\gamma }\left({\vec{\mu }}_{\alpha \gamma }{F}_{i,\gamma \beta }-F_{i,\alpha \gamma }{\vec{\mu }}_{\gamma \beta }\right)+\left(e\sum _{j=1}^N{Z}_j\frac{\partial {\vec{R}}_j}{\partial q_i}\right){\delta }_{\alpha \beta },$$

$$\frac{\partial }{\partial q_i}\vec{\mathbf{\mu }}=\left[\vec{\mathbf{\mu }},{\mathbf{F}}_i\right]+\left(e\sum _{j=1}^N{Z}_j\frac{\partial {\vec{R}}_j}{\partial q_i}\right)1,$$

$$\frac{\partial }{\partial q_j}{\vec{\mu }}_{21}\approx -\sum _{\gamma }{F}_{i,2\gamma }{\vec{\mu }}_{\gamma 1},$$

$${\psi }_t(q)=N_0\exp \left\{-(q-q_t{)}^{\top }\cdot A_t\cdot (q-q_t)+\frac{i}{\hslash }\left[p_t^{\top }\cdot (q-q_t)+{\gamma }_t\right]\right\}$$

$$V_{\text{eff}}(q,t)=V{|}_{q_t}+({\text{grad}}_{q}V{|}_{q_t}{)}^{\top }\cdot (q-q_t)+\frac{1}{2}(q-q_t{)}^{\top }\cdot {\text{Hess}}_{q}V{|}_{q_t}\cdot (q-q_t),$$

$${\dot{q}}_t$$

$${\dot{p}}_t$$

$${\dot{A}}_t$$

$${\dot{\gamma }}_t$$

$${\varphi }_0(q)=P\left(q-q_0\right){\psi }_0(q)=P\left(\frac{\hslash }{i}\frac{\partial }{\partial p_0}\right){\psi }_0(q).$$

$${\varphi }_t(q)=P\left(\frac{\hslash }{i}\frac{\partial }{\partial p_0}\right){\psi }_t(q),$$

$${\varphi }_t(q)=\left[\mu (q_0)+b_t^{\top }\cdot (q-q_t)\right]{\psi }_t(q).$$

$$b_t=\left(-2i\hslash A_t\cdot M_{qp,t}+M_{pp,t}\right)\cdot b_0,$$

Full text

{tocentry}

Vibrationally resolved electronic spectroscopy provides a valuable insight into the structure and dynamics of polyatomic molecules.1 Indeed, light-induced molecular dynamics is recognized as one of the key areas of research in physical chemistry, not only for fundamental understanding of Nature, but also for various applications, from solar cells to photodynamic therapy.2 The development of theoretical methods for simulating and understanding optical spectra is, therefore, of great importance.

The time-dependent approach to spectroscopy3 evaluates the vibronic spectrum as the Fourier transform of the nuclear wavepacket autocorrelation function, and, in contrast to the commonly used time-independent Franck-Condon approach,4, 5 can easily account for effects beyond the Born-Oppenheimer and global harmonic approximations. In the time-dependent approach, one must first perform exact or approximate molecular quantum dynamics. While the exact quantum dynamics typically requires a global potential energy surface6, 7 and thus scales exponentially with dimensions, semiclassical methods, such as the initial value representation,8 thawed Gaussian approximation (TGA),9 frozen Gaussian approximation,10 and Herman-Kluk propagator11, 12 require only local information and are suitable for on-the-fly implementation. The idea of using multiple frozen Gaussians as a basis for describing the full wavepacket has inspired a number of quantum13, 14, 15 and semiclassical 16, 17, 18, 19, 20, 21 “first-principles” approaches, which allow a combination with an on-the-fly ab initio (OTF-AI) evaluation of the electronic structure. Apart from a few examples,22 thawed Gaussians have been largely marginalized since they can neither describe wavepacket splitting nor very anharmonic dynamics. Yet, vibrationally resolved electronic spectra are mostly determined by short-time dynamics, during which both wavepacket splitting and anharmonic effects are less important. Indeed, a recent implementation of an on-the-fly ab initio thawed Gaussian approximation (OTF-AI-TGA) reproduced successfully vibrational structure of electronic absorption, emission, and photoelectron spectra, even in rather anharmonic and floppy systems such as ammonia. 23, 24

Here, an extension of the OTF-AI-TGA beyond the Condon approximation is presented by employing the Herzberg-Teller (HT) approximation,25 which, in contrast to the Condon approximation, 4 includes a linear dependence of the transition dipole moment on nuclear coordinates. We employ the OTF-AI implementation of the extended thawed Gaussian approximation (ETGA) 26 to evaluate the absorption spectra of phenyl radical, an anharmonic system allegedly exhibiting a significant Herzberg-Teller contribution, 27 and benzene, a textbook example of a symmetry-forbidden (i.e., electronically forbidden) transition.28

At zero temperature, within the electric-dipole approximation, first-order time-dependent perturbation theory, and rotating-wave approximation, the absorption cross section of a molecule with two electronic states that are not nonadiabatically coupled can be expressed as the Fourier transform

[TABLE]

of the dipole time autocorrelation function

[TABLE]

where $\hat{H}_{1}$ and $\hat{H}_{2}$ are nuclear Hamiltonian operators of the ground and excited electronic states, $\hat{\mu}$ is the matrix element $\hat{\vec{\mu}}$ ($=\hat{\vec{\mu}}_{21}$)111The subscript 21 is removed for simplicity since we will almost exclusively consider this matrix element. of the electric transition dipole moment operator projected along the polarization unit vector $\vec{\epsilon}$, i.e., $\hat{\mu}:=\hat{\vec{\mu}}\cdot\vec{\epsilon}$, and $|1,g\rangle$ is the ground vibrational state of the ground electronic state with zero-point energy $E_{1,g}$. Thus, the vibronic spectrum can be evaluated by propagating the initial wavepacket $|\phi(0)\rangle=\hat{\mu}|1,g\rangle$ on the excited state potential energy surface.

To compare with an experiment in gas phase, where the molecules are isotropically distributed, one must average the computed spectrum over all molecular orientations. Yet, due to the isotropy of space, a brute-force  numerical averaging is avoided; averaging is required only over three orientations of the molecule.29 To show this, consider the spectrum and autocorrelation function as $3\times 3$ tensors $\sigma_{\vec{\mu}\vec{\mu}}(\omega)$ and $C_{\vec{\mu}\vec{\mu}}(t)$, the latter related to $C_{\mu\mu}(t)$ by $C_{\mu\mu}(t)=\vec{\epsilon}^{\,\top}\cdot C_{\vec{\mu}\vec{\mu}}(t)\cdot\vec{\epsilon}\,.$ A simple analytical calculation30 shows that the average over all orientations of the polarization vector $\vec{\epsilon}$ is $\overline{C_{\mu\mu}(t)}=(1/3)\text{Tr}\left[C_{\vec{\mu}\vec{\mu}}(t)\right]$ and the corresponding absorption cross section is

[TABLE]

Thus, the average is easily evaluated, e.g., by averaging over only three arbitrary orthogonal molecular orientations with respect to the fixed polarization vector $\vec{\epsilon}$ or by fixing the molecular orientation and averaging over only three arbitrary orthogonal polarization vectors.

While the result (3) holds for arbitrary coordinate dependence of the transition dipole $\vec{\mu}(q)$, two approximations are frequently used. Within the most common Condon approximation, the transition dipole is considered constant: $\vec{\mu}(q)\approx\vec{\mu}(q_{0})$ and the general result (3) reduces to the textbook recipe for the averaged Franck-Condon (FC) spectrum:

[TABLE]

(“Divide by $3$ the spectrum for the molecular dipole aligned with the field.”) In the more accurate Herzberg-Teller approximation, the transition dipole moment becomes a linear function of nuclear coordinates:

[TABLE]

By explicitly differentiating the matrix element $\vec{\mu}_{\alpha\beta}$ of the molecular dipole between electronic states $\alpha$ and $\beta$ (for a moment we reintroduce the subscripts), one can show30 that

[TABLE]

or, in a more compact matrix notation:

[TABLE]

where $N$ is the number of atoms, $Z_{j}$ the atomic number and $\vec{R}_{j}$ the coordinates of the $j$th atom, and $F_{i,\alpha\beta}:=\langle\alpha|(\partial\beta/\partial q_{i})\rangle$ is the $i$th component of the nonadiabatic coupling vector between states $\alpha$ and $\beta$. For the transition dipole moment, $\alpha\neq\beta$ and the term proportional to $\delta_{\alpha\beta}$ in Eq. (6) [to $\mathbf{1}$ in Eq. (7)] vanishes, which shows that the Herzberg-Teller dependence originates from a combination of nonzero nonadiabatic and dipole couplings of states $\alpha$ and $\beta$ to an intermediate state $\gamma$. Moreover, since one may usually neglect nonadiabatic couplings between the ground and excited electronic states at the ground state optimized geometry, only the second term of the commutator survives:

[TABLE]

showing that a nonvanishing gradient of the transition dipole between the ground state $1$ and excited state $2$ requires an intermediate, “bright” ($\vec{\mu}_{\gamma 1}\neq 0$) excited state $\gamma$ that is vibronically coupled ($F_{i,2\gamma}\neq 0$) to the excited state $2$ . Such interpretation reveals the deep connection between the Herzberg-Teller approximation and the concepts of “vibronic coupling” and “intensity-borrowing.”31 Finally, although Eq. (7) suggests a way to evaluate $\partial_{q}\vec{\mu}_{\alpha\beta}$ from $\vec{\bm{\mu}}$ and $\mathbf{F}_{i}$, it is usually easier to evaluate the gradient by finite difference.

To find the autocorrelation function $C_{\mu\mu}(t)$, one must propagate the wavepacket. Heller’s TGA9, 3, 26 relies on the fact that a Gaussian wavepacket evolved in at most quadratic potential remains a Gaussian. Within this approximation, the anharmonicity of the potential is taken into account partially by propagating the Gaussian wavepacket

[TABLE]

in the time-dependent effective potential given by the local harmonic approximation of the full potential

[TABLE]

where $V|_{q_{t}}$, $\text{grad}_{q}V|_{q_{t}}$, and $\text{Hess}_{q}V|_{q_{t}}$ denote the potential energy, gradient, and Hessian evaluated at the center of the Gaussian, $N_{0}$ is the initial normalization constant, $(q_{t},p_{t})$ are the phase-space coordinates of the center of the Gaussian wavepacket at time $t$, $A_{t}$ is the complex symmetric width matrix, and $\gamma_{t}$ is a complex number; its real part gives an overall phase factor and its imaginary part ensures the normalization of the Gaussian wavepacket at all times. Parameters of the Gaussian follow Heller’s equations of motion9, 3, 26

[TABLE]

where $m$ is the diagonal mass matrix and $L_{t}$ the Lagrangian.

The extended TGA used in this work considers a more general form of the initial wavepacket, namely a Gaussian wavepacket (9) multiplied by a polynomial $P\left(q-q_{0}\right)$ in nuclear coordinates, which, at time zero, can be written as a polynomial in the derivatives with respect to $p_{0}$:

[TABLE]

This observation leads to a simple recipe for propagating the extended TGA wavepacket within the local harmonic approximation;26 namely, the wavepacket retains this form at all times:

[TABLE]

where $\psi_{t}(q)$ is the original TGA wavepacket (9). As for the initial Herzberg-Teller wavepacket, $P\left(q-q_{0}\right)=\mu(q_{0})+b_{0}^{\top}\cdot\left(q-q_{0}\right)$, where $b_{0}=$grad${}_{q}\mu|_{q_{0}}$, and the semiclassical propagation yields

[TABLE]

The four parameters of the Gaussian $\psi_{t}(q)$ are propagated with the usual TGA equations (11)-(14) and only one additional parameter, $b_{t}$, must be evaluated:

[TABLE]

where $M_{qp,t}=\partial q_{t}/\partial p_{0}$ and $M_{pp,t}=\partial p_{t}/\partial p_{0}$ are the elements of the stability matrix, which is already needed for propagating $A_{t}$ and $\gamma_{t}$, implying that the additional evaluation of $b_{t}$ comes at almost no additional cost. Remarkably, the orientational averaging of ETGA spectra is even simpler than what could be expected from the general simplification we mentioned in the beginning: within ETGA, the averaging requires only a single trajectory (instead of three) because the transition dipole moment does not affect the propagation of $q_{t}$, $p_{t}$, $A_{t}$, $\gamma_{t}$.

As a consequence of this property, the ETGA can be combined with an on-the-fly ab initio scheme at the same cost as the original TGA for spectra within Condon approximation.23, 24 The ab initio calculations are typically performed in Cartesian coordinates, and therefore the ab initio gradients and Hessians needed in Eq. (10) must be transformed23, 24 to the coordinate system $q$ that fits into our framework—the vibrational normal modes.

As the extended TGA is exact in a globally harmonic potential, it is useful to compare the on-the-fly approach with two common approximations of the excited-state potential energy surface: the vertical harmonic (VH) and adiabatic harmonic (AH) approximations,32 in which the excited-state potential is expanded to the second order about the ground and excited state optimized geometries, respectively (see Refs. 23  and 24 for details). We use density functional theory for the ground and time-dependent density functional theory for the excited state ab initio calculations; B3LYP/SNSD for phenyl radical and B3LYP/6-31+G(d,p) for benzene (see Supporting Information for details and validation by comparison with a higher level ab initio method).

According to Barone and coworkers,27, 33 the calculation of the absorption spectrum corresponding to $\tilde{\text{A}}^{2}\text{B}_{1}\leftarrow\tilde{\text{X}}^{2}\text{A}_{1}$ electronic transition of phenyl radical depends on the dimensionality of the simulation model and on the inclusion of the Herzberg-Teller contribution, anharmonicity effects, and mode-mixing (Duschinsky effect).34 Our model includes all these effects and provides the means to evaluate their importance.

To assess the influence of anharmonicity, the experimental spectrum is compared with the spectra simulated using the global harmonic approaches (Fig. 1, top). While the vertical harmonic approach only captures the overall envelope of the experimental spectrum, but fails to capture any details, the adiabatic harmonic model reproduces all main features of the spectrum. This is in contrast with the common expectation that the vertical harmonic approach should be more accurate,32, 35 as it describes better the Franck-Condon region of the excited state potential. Indeed, the emission spectra of oligothiophenes23 and both absorption and photoelectron spectra of ammonia24 are much better described with the VH than the AH approach. In phenyl radical, the failure of the VH approach lies in the incorrect frequencies and displacements of the two most displaced modes $18$ and $24$ (Fig. S4 of Supporting Information), resulting in the missing mode effect:36, 37 when the spectrum is not well resolved, it may contain a single progression whose spacing does not correspond to any of the vibrational frequencies of the system. Unlike the global harmonic approaches, the on-the-fly method overcomes the problem of guessing which excited state Hessian should be used (i.e., vertical or adiabatic) and reproduces the experimental spectrum rather well; moreover, with minimum human input.

To assess the validity of the Condon approximation, the FC and Franck-Condon Herzberg-Teller (FCHT) spectra of phenyl radical are compared at the bottom of Fig. 1. The FC and FCHT spectra are very similar since the absorption spectrum of phenyl radical is mostly determined by the symmetry-allowed Franck-Condon transition, while the Herzberg-Teller contribution only broadens the peaks slightly. We find that, in phenyl radical, including anharmonicity effects with the OTF-AI scheme is more important than including the HT contribution with the extension of the TGA.

Both the ground and excited state geometries of benzene belong to the $D_{6h}$ point group. Group theory predicts that the $\tilde{\text{A}}^{1}\text{B}_{2\text{u}}\leftarrow\tilde{\text{X}}^{1}\text{A}_{1\text{g}}$ electronic transition is symmetry-forbidden. Yet it is vibronically allowed, since the nonzero elements of the gradient of the transition dipole moment give rise to the vibronic spectrum.28 The gradient of the transition dipole moment originates mostly from nonadiabatic couplings between the $\text{B}_{2\text{u}}$ state and the bright $\text{E}_{1\text{u}}$ state.39 The spectrum contains one strong progression, assigned to one of the totally symmetric modes, as well as a number of hot bands. At present we do not attempt to simulate the hot bands, which require a finite-temperature treatment;40 our goal is computing the correct absorption cross sections of the main progression. In addition, we do not treat the spectral features arising from the second order vibronic coupling, i.e., a number of small intensity peaks which cannot be described within the first-order Herzberg-Teller approximation.41

Interestingly, the adiabatic harmonic approach again reproduces the experimental spectrum at least qualitatively, unlike the vertical harmonic approach which results in a number of peaks not observed in the experiment (Fig. 2, top). In the VH model, modes 25, 29 and 30 are more distorted since their frequencies are significantly lower than the corresponding frequencies in the AH model (see Table S6 and Fig. S5 of Supporting Information). Again, due to a partial treatment of anharmonicity, the OTF-AI-ETGA spectrum shows significant improvement over both global harmonic methods: While the relative intensities of AH peaks have errors of $20$ to $50\%$ and the VH model fails completely, the relative intensities of the OTF-AI-ETGA peaks lie within $5\%$ of experiment. Influence of anharmonicity on the spectrum is investigated using the autocorrelation functions in Supporting Information (see Figs. S6–S8). To further explore whether it is the error in the phase or magnitude of the autocorrelation function that affects the spectra more, we construct two hybrid, nonphysical autocorrelation functions. The first, denoted $|\text{AH}|~{}\exp(i\,\text{OTF})$, combines the magnitude of the adiabatic harmonic autocorrelation function with the phase from the on-the-fly correlation function, while the second, denoted $|\text{OTF}|~{}\exp(i\,\text{AH})$, combines the magnitude of the on-the-fly autocorrelation function with the phase from the adiabatic harmonic autocorrelation function. The spectra in Fig. 3 simulated using these two hybrid autocorrelations show clearly that it is the error in the phase of the AH model that most corrupts the intensities of the peaks. The difference between the two phases is closely related to the propagation of the stability matrix, and, consequently, to the Hessians of the excited electronic state potential.

In contrast to almost perfect description of the intensities by the OTF-AI-ETGA, neither the AH nor the OTF-AI-ETGA reproduces correctly the experimental spacing of the peaks in the main progression. This error, however, can be assigned to the electronic structure method, implying that the use of other density functional or wavefunction-based methods could give more accurate curvature of the potential, and hence a better spacing.

Our main result is contained in the bottom panel of Fig. 2, which compares the Franck-Condon spectrum (based on the Condon approximation), which can be evaluated with the original TGA, and the FCHT spectrum (based on the Herzberg-Teller approximation), which requires the extended TGA. While the FCHT agrees very well with experiment, the Franck-Condon spectrum is zero, since the transition dipole moment at the ground-state equilibrium is zero, which is the precise meaning of an electronically forbidden transition.

Incidentally, one can imagine simulating the FC spectrum by the commonly adopted procedure in which the transition dipole moment is set to unity and the spectrum is rescaled at the end. Although reasonable for electronically allowed transitions, here this “blind” procedure makes no sense since $\mu_{\text{FC}}=0$. This “FC” approach is compared to the FCHT result (see Fig. S9 of Supporting Information), and, somewhat surprisingly, the spectra are very similar, the main difference being a constant shift corresponding exactly to the excited state frequency of the degenerate inducing modes 27 and 28 (see Tables S5–S6 of Supporting Information). However, such similarity between the “FC” and the FCHT approach is not general. In addition, the “FC” approach is not capable of reproducing the absolute magnitudes of absorption cross sections, while the full FCHT approach provides a good estimate of the absolute absorption cross sections, as shown in Fig. S10 of Supporting Information.

To conclude, we presented an extension to the on-the-fly ab initio TGA and employed it to evaluate absorption spectra of phenyl radical and benzene within the Herzberg-Teller approximation. Considerable improvement compared to the usual global harmonic approaches was achieved and a further insight into the origins of the spectral features was given. The results obtained for the absorption spectrum of phenyl radical imply that including the anharmonicity effects is more important than the Herzberg-Teller contribution to the spectrum. Although thawed Gaussian approximation is often described in the context of calculating low-resolution spectra, here we reported the evaluation of a rather high-resolution absorption spectrum of benzene with high accuracy. The improvement was especially pronounced in the intensities of the peaks, due to a partial inclusion of the anharmonicity of the excited state potential. Thus, the OTF-AI-ETGA can be used not only to reproduce spectra, but more importantly, to evaluate the importance of different effects by going beyond the commonly used Condon and global harmonic approximations.

{suppinfo}

Computational details, ground and excited state optimized geometries, transition dipole moments, frequencies, and displacements, validation of the Born-Oppenheimer approximation and of the choice of the electronic structure method, conservation of the phase space volume and symplectic structure, details of spectra calculations, models describing the failure of the vertical harmonic method for absorption spectra, discussion of anharmonicity effects in benzene, “Franck-Condon” (“FC”) spectrum, and absolute absorption cross sections of benzene.

{acknowledgement}

The authors acknowledge the financial support from the Swiss National Science Foundation through the NCCR MUST (Molecular Ultrafast Science and Technology) Network, from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 683069 – MOLEQULE), and from the EPFL.