Photo-dissociation dynamics is simulated for vibrationally pre-excited pyrrole molecules using an *ab initio* multiple cloning approach. Total kinetic energy release (TKER) spectra and dissociation times are calculated. It is found that pre-excitation of N–H bond vibrations facilitates fast direct dissociation, which results in a significant increase in the high-energy wing of TKER spectra. The results are in very good agreement with the recent vibrationally mediated photo-dissociation experiment, where the TKER spectrum was measured for pyrrole molecules excited by a combination of IR and UV laser pulses. Calculations for other vibrational modes show that this effect is specific for N–H bond vibrations: Pre-excitation of other modes does not result in any significant changes in TKER spectra.

## I. INTRODUCTION

Ultrafast excited state dynamics is key to understanding many important processes in chemistry and biochemistry, such as light harvesting in plants^{1} and visual reception.^{2} It also plays an important role in the mechanisms of photoprotection,^{3} i.e., the processes that help biological molecules to minimize damage inflicted by UV light such as dissociation of the most important bonds.

Total kinetic energy release (TKER) spectra and velocity map images (VMIs) are efficient experimental tools that can provide valuable information about ultrafast processes following photo-absorption. In the TKER/VMI experiments, kinetic energy and velocity distribution of molecular fragments are measured. Recently, the VMI/TKER methodology technique has been modified and developed further in the vibrationally mediated photo-dissociation (VMP)^{4} method, where, prior to UV photoexcitation, the molecule is excited vibrationally. It has been demonstrated that IR excitation of certain vibrational modes of pyrrole prior to UV excitation changes the shape of its TKER spectrum.^{5} As not all vibrational modes have the same effect on the photo-dissociation, the VMP technique can be used to study mode selective photochemistry.

Heteroaromatic molecules, such as pyrrole, are the building blocks of many biomolecules, and their photo-dissociation has been extensively investigated by the above-mentioned experimental techniques. Photo-dissociation of pyrrole and its derivatives has been studied by TKER spectroscopy and time resolved VMI, yielding useful information about the mechanisms of dissociation and reaction rates. It has been shown that photo-dissociation mostly yields high kinetic energy hydrogen, which is scattered forward in the direction of the chemical bond. The energies of hydrogen, the time scale of its formation, and isotope effects have been measured.^{6–11}

Theoretically, the electronic structure of pyrrole is well understood since the seminal work of Sobolewski and Domcke^{12} where the role of πσ^{*} states in photo-dissociation was recognized. In a number of our previous studies, we applied our *Ab Initio* Multiple Cloning (AIMC)^{13,14} method to simulate the process of the photo-dissociation for pyrrole^{15–17} and other heterocyclic molecules.^{18,19} We demonstrated that the AIMC approach is capable of reproducing the main features of the experimental TKER spectra and VMI along with isotope effects providing a valuable insight into ultrafast non-adiabatic dynamics following the photo-excitation. In this paper, we apply the AIMC method to simulate photo-dissociation of vibrationally pre-excited pyrrole and compare our results with recent VMP experiment.^{5}

Electron-vibrational dynamics in molecules has been studied by many quantum techniques, which consider both electrons and nuclei quantum mechanically, such as *ab initio* multiple spawning^{20} and variational Multiconfigurational Gaussians (vMCGs).^{21} Recent examples also include on-the-fly semiclassical simulations of IR vibrational spectra^{22,23} and modeling of vibrationally resolved electronic spectra using the generating function approach^{24,25} and on-the-fly *ab initio* extended thawed Gaussian approximation.^{26,27} To the best of our knowledge, however, the process of vibrationally mediated photo-dissociation has never been simulated before.

In Sec. II, we briefly describe our AIMC approach. Section III discusses the initial conditions for the non-adiabatic dynamics with pre-excitation of particular vibrational modes. Section IV contains the computational details of our simulations. In Sec. V, we present and discuss the results. Conclusions are given in Sec. VI.

## II. AIMC METHOD

The *ab initio* multiple cloning (AIMC) approach belongs to the group of direct dynamics methods that run quantum non-adiabatic molecular dynamics using electronic structure data provided by “on-the-fly” calculations. The idea of AIMC is to combine the best features of Multiconfigurational Ehrenfest (MCE)^{28,29} and *Ab Initio* Multiple Spawning (AIMS)^{20} methods. Similar to MCE, the total wave-function |*Ψ*(**R**, **r**, *t*)⟩ is represented in a trajectory guided basis |*ψ*_{n}(**R**, **r**, *t*)⟩,

Each basis function consists of nuclear and electronic parts,

The nuclear part |*χ*_{n}(**R**, *t*)⟩ is a Gaussian coherent state,

where $R\xafn$ and $P\xafn$ are the coordinate and momentum of the center of *n*th Gaussian and the parameter *α* is the width of Gaussian functions. We use frozen Gaussians, where the width α is kept constant. Following Ref. 30, the value of α is taken as 4.7, 22.7, and 19.0 Bohr^{−2} for hydrogen, carbon, and nitrogen atoms, respectively. The electronic part of each configuration is represented in a basis of adiabatic states |*ϕ*_{I}(**r**; **R**)⟩ or electronic states in the center of Gaussian $\varphi Ir;R\xafn(t)$ in the Time Dependent Diabatic Basis (TDDB) version of the method.^{14,31} For sufficiently small molecules, such as pyrrole, the electronic eigenfunctions do not exhibit abrupt changes with the motion of the nuclei. In this case, the dynamics in both adiabatic and TDDB versions is described by the same set of equations presented below. See Ref. 13 for details.

The ansatz (2) is similar to the “single-set ansatz” used in the vMCG^{21} method. The difference is that in our approach the trajectories of different Gaussians are not coupled, with each of them guided by its own Ehrenfest force,

where *V*_{I}(**R**) is the *I*th potential energy surface and **C**_{IJ} is a non-adiabatic coupling vector **C**_{IJ} = ⟨*ϕ*_{I}|∇_{R}*ϕ*_{J}⟩. The Ehrenfest amplitudes *a*_{I} are propagated together with $R\xafn$ and $P\xafn$ as

where

Finally, coherent states (3) conventionally include phase *γ*_{n}(*t*) that is guided by the following equation:

Time dependence of the wave-function (1) is given by that of the trajectories $R\xafnt$ and $P\xafnt$, amplitudes $aInt$ in Ehrenfest configurations (2), and their amplitudes *c*_{n}(*t*) in the wave-function (1), which reflect the interaction between trajectories. An important advantage of the AIMC method is that the trajectories guiding basis functions are not coupled and can be calculated independently using Eqs. (4)–(7). Then, substituting (1) into the time dependent Schrödinger equation, we get the following equation for amplitudes *c*_{n}(*t*):

where *H*_{mn} are matrix elements of the Hamiltonian,

While the first term in (9) can be calculated analytically, the approximation is required for other terms. We use the first- and zeroth-order bra-ket averaged Taylor (BAT) expansion,^{13} the accuracy of which was recently demonstrated,^{32}

The principal advantage of BAT is that it utilizes only electronic structure data at trajectory points, which are needed anyway for basis propagation. As a result, propagation of amplitudes *c*_{n}(*t*) using BAT approximation can be performed in post-processing using saved trajectory data and does not require any additional electronic structure calculations, which are the most expensive part of the dynamics. Note that including the first-order term for potential energy matrix elements is extremely important, as it provides the imaginary part of the prefactor. More details can be found in Refs. 14 and 15.

Equations (4)–(11) form a complete set for calculating time evolution of the wave-function in the MCE approach. What distinguishes AIMC dynamics from MCE is cloning. Cloning is applied when several uncoupled electronic states are significantly populated: In this case, the average Ehrenfest trajectory becomes unphysical and should be branched to account for the wave-function bifurcation. In this case, a single Ehrenfest configuration is replaced by two new configurations: One of which has nonzero amplitudes for only one electronic state, while the second one contains contributions of all other electronic states,

and

The amplitudes of the two new configurations are adjusted in such a way that their total contribution to the whole wave-function (1) remains the same as the contribution of the original configuration,

Thus, cloning does not change the wave-function while providing additional flexibility for the ensuing dynamics.

Two different sets of criteria to initiate cloning have been proposed.^{13,33} In this work, as in our previous simulations for pyrrole^{15–17} and other heterocyclic molecules,^{18,19} we apply cloning when (1) the magnitude of breaking force

exceeds a threshold for at least one electronic state and, at the same time, (2) the magnitude of coupling between this state and all other states is below the second threshold. Thus, basis functions are usually cloned just after passing a conical intersection if, as a result of non-adiabatic population transfer, two electronic states with different potential energy gradients get significant Ehrenfest amplitudes.

AIMC is a fully quantum technique in a sense that it can be converged to the exact result, but in practice such convergence is not easy to achieve. We have developed and tested a number of tricks, which help us to approach convergence.^{19} As a first step, we expand the initial wave-function as a large sum of Gaussian “bits,” which are in turn represented in the basis of Gaussian Ehrenfest configurations (2). Then, due to the linearity of Schrödinger equation, instead of propagating the whole wave-function using a large basis, we propagate a large number of “bits” using their own small basis sets for each one of them. The minimal basis for such propagation would be just one trajectory guided Gaussian Ehrenfest configuration per the “bit” with clones originated from it. In Sec. III, we discuss the choice initial conditions for the “bits” used in this paper. Recently, we analyzed^{34} in great detail the convergence of AIMC and compared its performance with that of simpler but approximate techniques such as Ehrenfest dynamics (EHR) and Surface Hoping (SH), where it was demonstrated that in many cases the performance of AIMC is superior to that of EHR and SH. The improvements come at extra cost, which however can be greatly reduced by the use of sampling methods developed in AIMC.

## III. INITIAL CONDITIONS

Generating appropriate initial conditions for quantum non-adiabatic dynamics in trajectory-guided basis is a very non-trivial task that requires a number of approximations. In this paper, we use the simplest and the most straightforward approach assuming the Frank–Condon transition in which the electronic excitation simply “lifts” the ground electronic state vibrational wave-function into the exited electronic state. Although we recently proposed a more sophisticated approach^{17} where the process of electronic excitation is simulated in the course of dynamics, we use the Frank–Condon approximation here in order to reduce computational cost and ensure compatibility with our previous studies.^{15,16} Thus, the effect of vibrational pre-excitation will simply be modeled by using an excited vibrational state of ground electronic state as the initial condition for dynamics, which starts into the excited electronic state.

A set of random initial positions and momenta is normally generated from Wigner^{35} quasiprobability distribution for a harmonic oscillator using molecule Hessian in the ground state minimum. For the zero vibrational state, the Wigner distribution has the form

For exited vibrational states, however, the Wigner function cannot be used to generate the initial position and momenta, as it takes negative values around the zero. In particular, for the first vibrational state,

The obvious alternative is the Husimi Q distribution,^{36} which is always positive and can be seen as a smoothed version of the Wigner function.^{37} For zero and first states of harmonic oscillator, it has the form

Unfortunately, the Husimi function does not produce the correct expectation values of observables. In particular, the average classical energy for trajectories generated from the Husimi distribution would be *ℏω* instead of $12\u210f\omega $ and 2*ℏω* instead of $32\u210f\omega $ for zero and first vibrational states, respectively. This would not be too important if we had a near-complete basis, as in the case of model systems: AIMC is a fully quantum method, and the trajectories are used only to guide the basis to sample the most important part of the phase space at each particular moment of time. For real molecules, however, the basis in on-the-fly AIMC dynamics is always very far from being complete, and the excess trajectory energy in this case can affect the accuracy of the results. Taking all this into consideration, in this work, we generate initial conditions using the Husimi distribution but scale coordinates and momenta in order to get correct the trajectory energy. Using the scaling factor of $1/2$ for the zero state and $3/2$ for the first vibrational state, we get the following distribution:

## IV. COMPUTATIONAL DETAILS

Similar to our previous studies,^{15–19} simulations were performed using the AIMS-MOLPRO^{38} computational package modified to incorporate Ehrenfest dynamics. The complete active space self-consistent field (CASSCF) method at the SA4-CAS(8,7)/cc-pVDZ level was used for electronic structure calculations.

All trajectories are started from the S_{1} exited state with random initial positions and momenta generated from the probability distribution (19). Vibrational pre-excitation was taken into account by using distribution *P*_{1} for one particular vibrational mode with distribution *P*_{0} used for other modes.

An ensemble of 500 random trajectories was run for molecules without pre-excitation and with pre-excitation of the most important vibrational mode, which corresponds to N–H bond vibrations (mode 24). For pre-excitation of modes 15–23, a smaller ensemble of 250 trajectories was used because of the high computational cost of the calculations. We did not run calculations with pre-excitation of modes 1–14 taking into consideration that pre-excitation of low-frequency vibrations is unlikely to cause any significant changes in the TKER spectrum due to low energy surplus.

Trajectories were propagated for 200 fs, or until the N–H bond exceeded 3.5 Å, which was defined as the point of dissociation. As in our previous simulations of pyrrole,^{15} the cloning thresholds were taken as 5 · 10^{−6} Hartree/Bohr and 2 · 10^{−3} Bohr^{−1} for the magnitude of breaking force (15) and the norm of the non-adiabatic coupling vector, respectively.

Each trajectory originally guides one Gaussian with initial amplitude *c*_{n} = 1; the basis then grows as a result of cloning. For each branch leading to the dissociation, we calculate the H-atom kinetic energy *E*_{n} at the end of the trajectory *t* = *t*_{end} when it is sufficiently far and is no longer interacting with the radical. As the radical is much heavier than the hydrogen atom, *E*_{n} represent TKER of the dissociation. Solving Eq. (8) gives amplitudes *c*_{n}(*t*), which determine the weight of each branch contribution into the TKER spectrum, which is defined as |*c*_{n}(*t*_{end})|^{2}*δ*(*E* − *E*_{n}). The spectrum is then averaged over the ensemble and smoothed by replacing delta-functions with Gaussian functions (*σ* = 1200 cm^{−1}).

As S_{0} → S_{1} transition in symmetry prohibited in pyrrole, vibrational pre-excitation increases the oscillator strength and, therefore, the share of photo-excited molecules. To take this into account when comparing spectra with and without vibrational pre-excitation, the TKER spectrum for each vibrational mode was scaled according to the ratio of the ensemble-averaged square of transition dipole at initial positions for this mode to that without pre-excitation. The scaling factors for ten highest frequency modes are given in Table I.

Mode . | Scaling . | Mode . | Scaling . |
---|---|---|---|

15 | 1.165 | 20 | 1.143 |

16 | 1.076 | 21 | 0.956 |

17 | 1.717 | 22 | 1.240 |

18 | 1.120 | 23 | 0.923 |

19 | 1.140 | 24 | 1.277 |

Mode . | Scaling . | Mode . | Scaling . |
---|---|---|---|

15 | 1.165 | 20 | 1.143 |

16 | 1.076 | 21 | 0.956 |

17 | 1.717 | 22 | 1.240 |

18 | 1.120 | 23 | 0.923 |

19 | 1.140 | 24 | 1.277 |

In dissociation kinetics calculation, we first determine the so-called raw dissociation times *t*_{r} for all trajectories, which are defined as the times when the N–H bond length exceeds 3.5 Å threshold. Then, these raw times are blurred with Gaussian functions $GIRFt\u2212tr=exp\u2212t\u2212tr2/2\sigma XC2$ in order to account for the finite temporal widths of pump and probe pulses. We take here *σ*_{XC} = 51 fs, which corresponds to ∼120 fs FWHM of Instrument Response Function (IRF) in experiment.^{9} To obtain dissociation time constants *τ*, these dependencies are then fitted to the kinetic model^{9} *AG*_{IRF}(*t*) * [(1 − exp(−*t*/*τ*))*θ*(*t*)], where *θ*(*t*) is the Heaviside unit step function. See Refs. 18 and 19 for more details.

## V. RESULTS AND DISCUSSION

Calculated TKER spectra for the photo-dissociation of pyrrole with and without vibrational pre-excitation of mode 24 associated with N–H bond vibrations are presented in Fig. 2 (left). One can see that vibrational pre-excitation significantly increases the dissociation yield in the high-energy part of spectrum, while the low-energy part remains practically unchanged. These results are in extremely good agreement with experimental data for IR + UV photo-dissociation of pyrrole at *λ*_{UV} = 243 nm and *ν*_{IR} = 3532 cm^{−1} from Ref. 5 reproduced in Fig. 2 (right): The difference is small for low energies, drops to near zero in the middle of the left slope of the line, and then peaks with the maximum around the end of right slope. A small negative difference for near-zero energies in computational spectra probably can be attributed to the low accuracy in this area of the spectrum due to an insufficient number of low-energy trajectories. While the calculated energies are on average about 1.5 times higher than experimental values, which manifest itself in both the position and the width of the peak, this discrepancy can be ascribed to the lack of dynamic electron correlation in the CASSCF approach. The comparison of CASSCF and MS-CASPT2 energies for pyrrole indicates^{15} that the use of more accurate MS-CASPT2 potential energy surfaces would reduce the energies greatly improving quantitative agreement with experimental results.

Figure 3 compares the shapes of TKER spectra and the dissociation yield for the pre-excitation of different vibrational modes. One can see that only pre-excitation of N–H bond vibrations produces a noticeable high-energy wing in the spectra. The pre-excitation of any other modes, although typically increasing the integral yield of dissociation (mostly due to larger share of photo-excited molecules), does not lead to any significant changes in the shape of spectrum.

Along with TKER spectra, we also study the effect of vibrational pre-excitation on dissociation kinetics. Figure 4 presents normalized dissociation counts (both raw and smoothed with IRF) for molecules with and without vibrational pre-excitation of the N–H bond. Similar to our previous calculation,^{18,19} raw data exhibit a short delay and a very steep increase in dissociation count immediately afterward due to trajectories that dissociate right away or just after one vibrational period. Then, the dissociation significantly slows down, as the remaining trajectories must first sample the potential energy surface to find a way out.

Vibrational pre-excitation leads to significantly faster dissociation. Fitting our results with the kinetic model^{9} gives pyrrole dissociation time constants of 29.5 fs for pre-exited molecules and 48.2 fs without pre-excitation. The latter constant is in good agreement with experimental results of 46 ± 22 fs for *λ*_{UV} = 238 nm pyrrole photo-dissociation.^{9} To the best of our knowledge, no experimental measurements of dissociation time constant for vibrationally pre-exited pyrrole have yet been performed.

More details of dissociation kinetics can be seen from non-normalized data plotted in Fig. 5. These data show that vibrational pre-excitation leads to a large number of additional dissociations in the first few femtoseconds of dynamics, with the maximum around 15 fs time. On the other hand, the rate of dissociation at later times remains practically the same. This is consistent with TKER spectrum data, as the trajectory that dissociates right away is mostly responsible for a high energy part of TKER spectra.

## VI. CONCLUSIONS

We simulate the process of photo-dissociation of pyrrole with pre-excitation of different vibrational modes and compare results with the TKER spectrum from IR + IV experiment.^{5} We use the AIMC computational approach that has proven its efficiency in our previous simulations for pyrrole and other heterocyclic molecules. We show that our simulations for pyrrole with pre-excitation of N–H bond vibrations reproduce qualitatively the difference between UV and IR + UV TKER spectra observed in the experiment. We demonstrate the selectivity of vibrational pre-excitation, as pre-excitation of no other mode leads to similar changes in the TKER spectrum. Our calculations of photo-dissociation kinetics predict that pre-excitation of N–H bond vibrations should result in a lower dissociation time due to a significant increase in the dissociation rate in the first moments after electronic excitation, while the rate at later times remains practically unchanged.

## ACKNOWLEDGMENTS

We acknowledge the support from EPSRC through Grant No. EP/P021123/1. This work was undertaken on ARC3 and ARC4, part of the High Performance Computing facilities at the University of Leeds, UK.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.