Simulation of the dynamics of vibrationally mediated photodissociation for deuterated pyrrole

The dynamics of photodissociation for vibrationally pre-excited deuterated pyrrole molecules is simulated using ab initio multiple cloning (AIMC) approach. Total kinetic energy release (TKER) spectra and dissociation times are calculated. The results for pyrrole and deuterated pyrrole molecules with and without vibrational pre-excitation are compared. Calculations show that, as expected, the kinetic energy of additional dissociation fragments is lower in deuterated pyrrole and mostly located in the upper-middle part of the TKER spectrum. However, despite lower energy of dissociative bond vibrations, pre-excitation of deuterated pyrrole leads to higher dissociation yield increase than in pyrrole and significantly shortens dissociation time.


Introduction
The study of ultrafast excited state dynamics is essential to understand many fundamental biochemical processes. Photosynthesis in plants, visual reception of cone cells in the retina and the photostability of biological molecules, such as DNA, all take place by photochemical reactions, utilising electronic excitation and nonadiabatic transitions. In photosynthetic plants and organisms, non-adiabatic transitions are present in protein pumps which move electrons through a protein chain to a target area. The photoprotection mechanism in DNA uses nonadiabatic transitions to prevent damage of DNA molecules from UVA and UVB light, which would otherwise result in mutations leading to cancer and other genetic deformations. Pyrrole is a key building block that is present in a large number of these biological molecules, including DNA and chlorophyll, and the dynamics of its photodissociation (along with electron impact dissociation) has been a subject of extensive experimental [1][2][3][4][5][6][7][8][9][10][11][12][13] and theoretical [14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31] studies. Photofragment translational spectroscopy [32] is among the most efficient experimental tools used to study the process of photodissociation. We previously demonstrated [23,33,34] that simulations using our ab initio multiple cloning (AIMC) approach [35,36] can reproduce the principal features of experimental total kinetic energy release (TKER) spectra and velocity map images (VMI) for the photodissociation of pyrrole and other heterocyclic molecules, including the isotopic effect for their deuterated isotopologues. The dynamics calculations add valuable novel insight into the dissociation mechanism providing details that cannot be seen directly in the experiment.
Recently we proposed [31] that, with appropriate initial conditions, AIMC approach can be used to simulate the process of vibrationally mediated photodissociation (VMP) [37][38][39], when the molecule is vibrationally preexited by IR laser pulse before electronic excitation by UL pulse. Our calculations for vibrationally pre-excited pyrrole using this approach produced TKER spectra that show good qualitative agreement with VMP experiment [12]. In this work, we report our simulations of photodissociation dynamics for vibrationally preexcited deuterated pyrrole molecules.
In section 2, we briefly outline AIMC computational approach and specify the initial conditions for dynamics with and without vibrational pre-excitation. Section 3 provides computational details of our simulations. In section 4, we present and discuss the results. Conclusions are given in section 5. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Theory
The process of photodissociation is simulated using AIMC quantum non-adiabatic molecular dynamics method, which was already described in details in a number of our previous works [23,31,35,36]. Below we just briefly outline the main features of AIMC approach.
The idea of AIMC is to represent the total electron-nuclear wave-function in a basis of Gaussian coherent states moving along Ehrenfest trajectories: where R and r are nuclear and electron coordinates respectively. The electronic part is represented either in the basis of R-dependent adiabatic states f f ñ º ñ r r R ; I I | ( ) | ( ) or, in the alternative version of the method [40], in time-dependent diabatic basis (TDDB) f f ñ º ñ t r R ; I I n | | (¯( )) that is changing with the motion of Gaussians. It was shown [35] that for small molecules without abrupt changes of electronic eigenstates, such as pyrrole, the dynamics for both adiabatic and TDDB ansatzes is described by exactly the same set of equation.
The equations of motion for Gaussians are solved simultaneously with equation for the evolutions of Ehrenfest amplitudes: where is H IJ n ( ) is the electronic Hamiltonian for nth trajectory: and V I and C IJ are potential energies and non-adiabating coupling matrix elements (NACME).
The motion of Gaussians is guided by force: where the first term is the usual Ehrenfest state-average force and the second term, often called the Hellmann-Feynman force, is associated with non-adiabatic electronic population transfer. The Gaussians are moving independently from each other, and the coupling between trajectories is reflected only in the evolution of amplitudes D t , n ( ) which is obtained by solve time-dependent Schrödinger equation in basis (1). Ehrenfest trajectories become unphysical when two or more uncoupled electronic states with significantly different gradients have large amplitudes. In order to address this issue, we apply cloning procedure: we replace one 'mixed' basis function by two new ones, each dominated by just one electronic state. Both new basis functions have the same Gaussian part and, due to the appropriate choice of new amplitudes, their combined contribution to the total wave-function (1) remains the same as that of the original basis function. The state is cloned out when its breaking force exceeds a threshold and, at the same time, the magnitude of coupling between this state and all other states is below the second threshold. Cloning results in the branching of trajectories, which reflects the bifurcation of wave-function on conical intersections.
In order to generate the initial conditions for the dynamics, we assume that ground electronic state vibrational wave-function is simply 'lifted' into the exited electronic state. So, vibrational frequencies and normal modes are calculated for ground state potential energy surface, which are used to generate a random set of initial positions and momenta using quasiprobability distribution for harmonic oscillator. The effect of vibrational pre-excitation is modelled by using an excited state vibrational wave-function for a particular mode. As Wigner distribution takes negative values for excited vibrational states, it cannot be used here and, similar to our previous simulations for vibrationally pre-excited pyrrole [31], we use Husimi distribution which was appropriately scaled to provide correct average energy for both ground and first excited vibrational states: One can see that distribution (5) gives correct classical energy for trajectories, which is important considering extremely small size of the basis used in on-the-fly AIMC quantum dynamics. For ground state, P x p , 0 ( ) coincides with Wigner distribution providing compatibility with our previous simulations of pyrrole photodissociation dynamics [23-25].

Computational details
As before [23-25, 31, 33-35], the simulations were run using AIMS-MOLPRO [41] computational package modified to incorporate Ehrenfest dynamics with electronic structure date provided by complete active space self-consistent field (CASSCF) calculations at SA4-CAS(8,7)/cc-pVDZ level. The CASSCF method typically overestimates the excitation energies due to the lack of dynamic electron correlation. As a result, calculated kinetic energies are on average about 1.5 times higher than experimental values, which leads to the shift of calculated TKER spectrum. Nevertheless, as we show before [23], the AIMC calculations using CASSCF still correctly reproduce the main features of experimental spectra.
The trajectories are started from first exited state with random initial positions and momenta generated using distribution P x p , 1 ( ) for pre-exited vibrational mode and P x p , 0 ( )for all other modes. Then, they were propagated for 200 fs, or until the N-D bond exceeded 3.5 Å, which was defined as the point of dissociation.
Each trajectory is initially guiding a basis of one Gaussian coherent state with width α taken [42] equal to 4.7, 22.7 and 19.0 Bohr −2 for hydrogen, carbon and nitrogen atoms respectively. The basis is growing at cloning points and the evolution of amplitudes D n reflects the interaction between branches. As before, we used values of -5 10 6 · atomic units and -2 10 3 · atomic units as cloning thresholds for the magnitudes of breaking force and non-adiabatic coupling respectively.
In TKER spectra calculations, we find kinetic energy at the end of each trajectory branch leading to dissociation. The weight of each branch contribution into TKER spectrum is defined as D n 2 | | in the final point. The spectrum is then averaged over an ensemble and smoothed using Gaussian functions with s = -1200 cm . 1 We use an ensemble of 500 trajectories for molecules without vibrational pre-excitation and with pre-excitation of dissociative bond vibrations. According to our previous experience, this number of trajectories is optimal for TKER spectrum calculation: using larger ensemble, although computationally expensive, does not significantly improve the accuracy of the spectrum. For pre-excitation of other vibrational modes, we use a smaller ensemble of 300 trajectories, which is nevertheless sufficient to provide more or less correct shape of the spectrum.
For all trajectories, we determine raw dissociation times t , r which are defined as the times when N-D bond length exceeds the dissociation threshold. To calculate dissociation kinetics, these times are then smoothed with Instrument Response Function (IRF) ) which reflects the temporal widths of pump and probe pulses in the experiment. For the comparability with our results for H-pyrrole, we take here the dissociation threshold of 3.5 Å and s = 51 XC fs, the same as in our previous simulations [31]. The dissociation time constants t is then calculated by fitting the dissociation kinetic to the model [9] where q t ( ) is Heaviside unit step function. More details can be found in [33,34]. Figure 1 compares calculated TKER spectra of pyrrole and deuterated pyrrole (referred below as H-pyrrole and D-pyrrole). As one would expect, the spectrum for D-pyrrole has more intensive low-energy wing and its maximum is shifted by ∼1000 cm −1 in the low-energy directions. This effect of deuteration on the TKER spectrum is similar to our previous results for imidazole and pyrazole [34], although significantly more pronounced in the case of pyrrole. Figure 2 presents TKER spectrum for D-pyrrole with pre-excitation of vibrational mode 20, which corresponds to N-D vibrations. For comparison, we also show the spectrum for H-pyrrole with mode 24 preexcitation (N-H vibrations) from our previous work [31]. In order to compare spectra with and without vibrational pre-excitation, we must consider that pre-excitation also increases the oscillator strength (and, as a result, the share of photo-excited molecules) for symmetry prohibited S 0 →S 1 transition. This is taken into account by scaling the spectra for pre-excited molecules by 1.277 for H-pyrrole and by 1.119 for D-pyrrole, Figure 2. Calculated TKER spectra for the photodissociation of pyrrole (top) and deuterated pyrrole (middle) with and without preexcitation of dissociative bond vibrations. The integral yield is normalized to 1 for both spectra without vibrational pre-excitation. The difference between spectra with and without vibrational pre-excitation (bottom).

Results and discussion
which is the ratio of the ensemble-averaged square of transition dipole momenta for molecules with and without pre-excitation calculated at the starting points of the trajectories.
As one can see both for H-pyrrole and D-pyrrole, vibrational pre-excitation results in significant increase in the dissociation yield in high-energy part of TKER spectrum while low-energy part does not change much. However, the additional yield for D-pyrrole is mostly located around the maximum and the right slope of the peak, while for H-pyrrole it also extended to the area of much higher energies where the spectrum is close to zero without pre-excitation. This can be explained by lower energy of dissociative bond vibrations in D-pyrrole. Nevertheless, despite lower vibrational energy, the integral increase of dissociation yield is surprisingly higher for D-pyrrole and reaches 55%, versus 51% for H-pyrrole.
The selectivity of vibrational pre-excitation for the photodissociation of D-pyrrole is demonstrated in figure 3, which compares TKER spectra for molecules with pre-excitation of vibrational modes 18-22, where mode 20 corresponds to N-D bond vibrations. Due to expensiveness of calculations, we have not run the dynamics for other vibrational modes. One can see that only pre-excitation of mode 20 results in the long intensive high-energy wing in the spectrum, while the pre-excitation of other modes does not cause significant changes in the spectrum shape. The same kind of selectivity was previously found [31] for H-pyrrole photodissociation, where only pre-excitation of mode 24 corresponding to N-H vibrations resulted in similar changes of TKER spectrum.
The dissociation kinetic is presented in figure 4, which compares normalized dissociation counts for H-pyrrole and D-pyrrole molecules with and without vibrational pre-excitation of the dissociative bond. From the insert showing raw dissociation counts, one can see that, similar to H-pyrrole, the photodissociation of D-pyrrole is a combination of fast and slow processes: some molecules dissociate directly and the others must first sample the potential energy surface to find a way. The pre-excitation of dissociative bond vibrations creates favourable conditions for fast direct dissociation significantly reducing the dissociation time. This effect is more pronounced for D-pyrrole, where the vibrational pre-excitation reduces dissociation time from 93.7 fs to 46.5 fs, versus 48.2 fs and 29.5 fs for H-pyrrole [31]. Considering lower surplus of vibrational energy in deuterated pyrrole, this behaviour looks unexpected. This computational prediction is waiting for experimental confirmation, as, to the best of our knowledge, the deuterated pyrrole has not been yet the subject of VMP study.

Conclusion
We have simulated the dynamics of vibrationally mediated photodissociation for deuterated pyrrole using AIMC computational approach and compare the results with our previous similar calculations for pyrrole. As in the case of pyrrole, calculations show that pre-excitation of N-D bond vibrations in deuterated pyrrole significantly facilitates the process of fast direct photodissociation producing more high-energy fragments. This reduces more that by half the dissociation time for deuterated pyrrole, significantly more than for pyrrole. The integral increase of the dissociation yield in also higher for deuterated pyrrole despite the lower frequency of dissociative bond vibrations, which manifest itself only in significantly less extended high-energy wing of its TKER spectrum. These two features are somewhat counterintuitive. Calculations for different vibrational modes demonstrate the selectivity of vibrationally mediated photodissociation, as only the pre-excitation of dissociative bond vibrations induces noticeable changes in the shape the TKER spectrum. The unexpectedly strong increase of dissociation yield in deuterated pyrrole along with strong decrease of its dissociation time, larger than in nondeuterated molecule, are unexpected. We hope that our numerical results may encourage future VMP experimental studies for deuterated pyrrole and other similar molecules.