Simulation of ultrafast photodynamics of pyrrole with a multiconfigurational Ehrenfest method

Abstract

We report the first results of ab initio multiconfigurational Ehrenfest simulations of pyrrole photodynamics. We note that, in addition to the two intersections of 1¹A₂ and 1¹B₁ states with the ground state 1¹A₁, which are known to be responsible for N–H bond fission, another intersection between the 1²A₂ and 1²B₁ states of the resulting molecular radical becomes important after the departure of the H atom. This intersection, which is effectively between the two lowest electronic states of the pyrrolyl radical, may play a significant role in explaining the branching ratio between the two states observed experimentally. The exchange of population between the two states of pyrrolyl occurs on a longer scale than that of N–H bond fission.

---

1. Introduction

Heteroaromatic molecules such as pyrrole, imidazole, and phenol have been the focus of experimental and theoretical attention due to the role their derivatives play as the chromophores of nucleobases and aromatic amino acids. The mechanisms of their photochemistry have been studied experimentally^1,2 and theoretically,^3–10 where it was suggested that the N–H/O–H bond fission was an important channel in the photodissociation dynamics. The role of the ¹πσ* states in this process has been emphasized^11–23 wherein a number of total kinetic energy release (TKER) spectra for various bond fission reactions were reported.

In particular, the photodissociation of pyrrole (C₄NH₅) has been investigated^24,25 and has proven to be an interesting case, prototypical of many heteroaromatic systems in which low energy heteroatom centred πσ* states occur in the singlet manifold. These states result from the close degeneracy of the 3s Rydberg orbital of the heteroatom with the repulsive σ* orbital of the heteroatom to substituent bond, which leads to strong mixing. Where the substituent is simply H, N–H bond fission is the primary decay mechanism for pyrrole molecules excited to this πσ* state. Experimental evidence for this dissociation has been reported²⁶ wherein H atoms were detected with a centre-of-mass frame energy of 18 kcal mol⁻¹ and a strong anisotropy, indicating prompt dissociation (within the molecular rotation period). Further evidence from ion-imaging¹ and hydrogen (Rydberg) atom photofragment translational spectroscopy (HRAPTS)²⁴ supported these initial findings. Additional details from the high resolution HRAPTS data showed several notable features, including direct (Herzberg–Teller allowed) excitation to the ¹A₂ electronic state, with the preservation of the H–T active vibrations in the product radical, apparently showing strong vibrational adiabaticity during the dissociation. Theoretical work by Sobolewski and Domcke^3,4 showed the existence of conical intersections in the decay channel, which are by now well known and opened the possibility of electronic population branching during the dissociation. Despite sufficient total energy to access both the ²A₂ and ²B₁ states of pyrrolyl, the interpretation of the HRAPTS data is that population was only observed in the lowest (X²A₂) electronic state of the radical. The peaks in the TKER spectrum were assigned to various vibrational states of this ground electronic state of the pyrrolyl radical.

Time resolved studies have also been performed^27,28 albeit in the region of the stronger¹ ππ* absorption, and place an upper bound on the timescale for N–H fission of ∼200 fs, although direct dissociation on the ¹πσ* state is expected to be much faster. Previously, quantum dynamics calculations^1,5,6 have been performed on a limited dimensionality model PES and have predicted rapid dissociation and significant population of the ²B₁ state of pyrrolyl, in contrast to the experimental data.

In this article we report the results of our simulations of the nonadiabatic dynamics of pyrrole with the help of the multiconfigurational Ehrenfest (MCE) method. In this approach, the nonadiabatic dynamics are represented by a trajectory which is “shared” between several electronic states. In the MCE^29–32 method, Ehrenfest trajectories are dressed with Gaussian wave packets which serve as a basis for exact quantum propagation of the nuclear wave packet. They carry and exchange quantum amplitudes, which makes MCE fully quantum in principle. However, as has been explained in our previous work, with time the Ehrenfest trajectories start running away from each other and the quality of the basis deteriorates. As a result, the quantum coupling between Gaussians decreases, although this can be corrected, in principle, by increasing the number of basis trajectories. Nevertheless even without full, long time quantum coupling between the amplitudes of trajectories, many features of the dynamics can be understood (at least qualitatively) on the basis of the analysis of an ensemble of almost independent Ehrenfest trajectories. This analysis provides a semi-classical approximation to the exact quantum dynamics of the nuclei, whilst still retaining a fully quantum description of electronically nonadiabatic processes. An important feature of the ab initio MCE method described in ref. 30 is that it is a direct dynamics which, similar to other ab initio techniques, calculates electronic potential energy surfaces “on the fly” and does not require any a priori information about the PESs. This is similar in spirit to other existing methods such as the ab initio multiple spawning^33–36 and ab initio variational Gaussians.^37,38 These alternative quantum direct dynamics methods differ from MCE and from each other only in the precise way the trajectories are guided.

The availability of “on-the-fly” ab initio Ehrenfest quantum dynamics allows a more complete analysis of the dissociation of pyrrole. By including all nuclear degrees of freedom in our numerical experiment, we do not bias the simulation by the choice of dimensions to a model, hence allowing any nuclear coordinate to influence the dynamics. The current simulation provides evidence that the population of the ²B₁ state is low and thus confirms the original assignment of the HRAPTS data. We notice that after partial departure of the H-atom, the two electronic states ²A₂ and ²B₁ are still coupled via an intersection between the potential energy surfaces of the pyrrolyl radical, which is nearly isolated. This intersection is reached by the pyrrolyl vibrational motion, and as a result the population of the higher ²B₁ state resulting from rapid nonadiabatic dynamics of the N–H bond fission leaks to the lower ²A₂ state via longer time dynamics at this final intersection. A possible role for this feature has been proposed in ref. 5 but a multidimensional simulation allows one to see it very naturally. Previously, such a conical intersection between two lowest electronic states of pyrrolyl has been found in ref. 39 and our current simulation suggests that it may play an important role in the nonadiabatic dynamics of pyrrole, which to the best of our knowledge has never been considered before.

2. Theory

Our multiconfigurational Ehrenfest approach to nonadiabatic dynamics is described in ref. 29–32. In particular see ref. 30 for ab initio MCE. In this article we only sketch briefly the main concept. The exact molecular wave function at time t can be expanded as a sum of N Ehrenfest wavepackets with amplitude D_k(t)The basis wavepackets |ψ_k(t)〉 are called “configurations” and the wave function (1) is therefore multiconfigurational. Each Ehrenfest configuration includes electronic and nuclear parts:The nuclear motion is described by a single trajectory dressed with a frozen Gaussian (FG) , which is a product of 3D frozen Gaussians for each atom j included in the quantum dynamics simulation. Often a frozen Gaussian is called a coherent state |z(t)〉 and labelled with a single complex vector

z_3D = (γ^1/2q_3D + iγ^−1/2ℏ⁻¹p_3D)/2^1/2 (3)

which includes both the real momentum and position of a particle. The parameter γ describes the width of a 3D Gaussian wave packet. The FG coherent state is “shared” between the two (or more) electronic states , with the amplitudes of n_e electronic states. A single configuration (2) is not flexible enough to describe accurately the quantum dynamics and it is well accepted that the standard Ehrenfest method is often unsuitable for accurate nonadiabatic dynamics. However the Ehrenfest configurations can be used as a basis function for a multiconfigurational approach and the wavefunction (1) can be converged to the exact result.^27,28

The evolution of the wave function is described by the motion of Ehrenfest configurations (z,a) and time dependence of their amplitudes D. The equations of motion for the Ehrenfest trajectories follow from the Hellmann–Feynman theorem.⁴⁰ They can be written in the form of Newton's second law, ṗ = F and q̇ = p/m, where

is the Ehrenfest force in the adiabatic representation. After the amplitudes are represented as a product of rapidly oscillating exponents of the classical action and the smooth pre-exponential factor, where , the equations for the amplitudes take the formwhere

Δ^(lj) = q̇δ^(lj)(q) (6)

and δ^(lj) = 〈ϕ^(l)|∇|ϕ^(j)〉, and δ^(jl) = 〈ϕ^(j)|∇|ϕ^(l)〉 = −δ^(lj) are the nonadiabatic coupling matrix elements (NACMEs) which couple the electronic states in the adiabatic representation. Conveniently, Ehrenfest trajectories are independent of each other.

The quantum coupled equations for the amplitudes D_k(t) in (1) are as follows:

The matrix elements of the Hamiltonian include those of kinetic energywhich are calculated analytically and potential energy 〈z_l|V^(jj)|z_k〉. For the latter we are using the linear approximationand a similar approximation for the NACMEswhich couple the adiabatic states V^(jj). An important feature of this approximation is that for quantum coupled equations it uses the information already available from Ehrenfest trajectories, so that quantum dynamics comes at almost no extra cost to ab initio Ehrenfest trajectories. Nonadiabatic coupling at conical intersections is singular, which cannot be captured by the simple linear approximation (10). However this singularity is important only if the conical intersection is between the two coupled trajectories. In multidimensional systems a swarm of trajectories never “hits” the singularity point. Only the trajectories which are very close in phase space are coupled and the chance of them passing on two different sides of the intersection is negligible. To verify the validity of approximation (10) we compared matrix elements (10) with those obtained by centroid approximation, which uses NACMEs in the middle between trajectories and found almost no difference for the trajectories close enough to be coupled.

The adiabatic potential energy surfaces V^(jj), their derivatives and the NACME vectors can easily be calculated using MOLPRO⁴¹ or any other electronic structure package.

3. Computational details

3.1 Initial sampling

The initial conditions for the Ehrenfest trajectories are based on projection of a ground state zero-point nuclear eigenfunction onto the S₁ surface, allowing excitation at all energies. Therefore, the simulated wavepackets contain a broader spread of total energies than is accessible in the experiments.²⁴ Initial conditions for the trajectories were generated from a vibrational Wigner distribution⁴² of the v = 0 vibrational state of the S₀ electronic stateBy performing the calculation of normal modes with associated frequencies and reduced masses, the phase-space distribution function for the product of the Wigner distributions was generated for each normal mode. The initial q and p are generated from this distribution (v = 0) using a Monte-Carlo sampling procedure. This sampling was also used by Martinez.³⁴ With sufficient statistics, MCE Ehrenfest configurations exchange their amplitude as described by the equations for the coefficients D in (1). However, if the number of trajectories is small, then they run away from each other quickly. After that point, quantum coupling (7) disappears and coefficients D become constant, the configurations are effectively uncoupled and MCE works as a semiclassical approximation to the exact quantum solution for nuclear dynamics, similar to the Frozen Gaussian approximation, which still gives a good qualitative (and sometimes even quantitative) description of the dynamics. To retain coupling between Ehrenfest configurations for longer and at the same time to reduce computational cost, we applied a strategy of propagating the wave function “bit by bit”, which has been used by us previously.^31,43 First, we ran several (namely 24) trial trajectories with randomly selected initial conditions (11). Only for 3 of them the N–H bond fission occurred quickly. For those 3 trajectories, which represent the dissociative “bits” of the wave function, we generated a bundle of 10 trajectories with very close initial conditions, which were used as initial basis. This procedure of propagating only important parts of the wave function “bit-by-bit” has been shown to improve greatly the quality of MCE propagation.³¹

3.2 The level of electronic structure theory

In this study, Dunning's cc-pVDZ (correlation-constant polarized valence double zeta) basis set⁴⁴ was employed on all atoms. In order to describe the low-lying πσ* states (1¹A₂ and 1¹B₁ states), one additional diffuse s function, one additional set of p functions and one additional set of d functions were added to the nitrogen atom, and one additional diffuse s function and one additional set of p functions were added to the dissociative hydrogen atom. The exponents of these orbitals were taken from the diffuse functions of Dunning's aug-cc-pVDZ basis set.⁴⁴ The necessary energy, gradient and NACME data for the S₀, S₁ and S₂ states were computed on the fly at the CASSCF level via interface with MOLPRO 2010.1.⁴¹ The complete active space consisted of seven orbitals (three ring π orbitals and two corresponding π* orbitals, one σ orbital and a corresponding σ* orbital) with eight electrons, and state averaging was performed over four states (SA4-CAS(8,7)-SCF). This level of theory is required to realistically represent the multiconfigurational nature of the electronic states as the N–H bond is extended and to allow analytic calculation of the non-adiabatic coupling between states, particularly near points of degeneracy. This level of theory generally reproduces the lowest three PESs for pyrrole, although the minimum in the S₁ surface due to the 3s Rydberg character of the ¹πσ* state is somewhat overestimated, leading to some artificial trapping of the wavepackets in this region. However, as in this study we will discuss the dynamics of the dissociative surfaces, the artificially deep well around the Franck–Condon region has little impact on the dynamics, once the dissociative part of the potential is accessed. Calculation of the electronic surfaces and conical intersection in pyrrolyl used identical methodology, save for the absence of the dissociative H atom and its single electron. Furthermore, only three states exist at low energies, leading to the SA3-CAS(7,7)-SCF calculation. For consistency, the same cc-pVDZ basis set, with additional basis functions applied to the N atom, was used in these calculations.

4. Results and discussion

Tables 1 and 2 show the orbitals which are taken into account in the active space of the CASSCF calculation, and the primary electron configurations of the lowest three electronic states at the equilibrium geometry (r_NH = 0.996 Å) and in the region H-atom dissociated (r_NH = 5 Å). The equilibrium geometry of pyrrole has C_2v symmetry. Fig. 1(a) illustrates the one-dimensional potential energy cuts of the low-lying electronic states as a function of the N–H stretching coordinate with a fixed molecular framework. As discussed in ref. 6 the CASSCF method underestimates the vertical excitation energies^45–49 so that the lines are plotted with an offset of 0.8 eV in order to compare with the CASSCF surfaces corrected by the MRCI method (plotted as dotted curves).⁶ The 1¹A₂ and 1¹B₁ curves calculated by the pure CASSCF method are in good agreement with those by the CASSCF/MRCI in the Franck–Condon region. Fig. 1(b) shows another 1D cut of the potential energy surfaces as a function of the N–H distance with a relaxed molecular skeleton. In the dissociative region the Ehrenfest trajectories of basis CSs tend to propagate in the vicinity of this minimum energy path region rather than the “frozen” PESs shown in Fig. 1(a). The states 1¹A₁, 1¹A₂, 1¹B₁, 2¹A₁ and 2¹B₁ are depicted by blue, red, green, orange and black lines. However, in our simulation, molecules no longer keep their symmetry and the states should be labelled as 1¹A, 2¹A, 3¹A, and 4¹A under C₁ symmetry, or with adiabatic labels: S₀, S₁, S₂, and S₃. Therefore Fig. 1 also indicates the S₀, S₁, S₂, and S₃ states, with dashed, solid, dotted and dash-dotted lines respectively. In all figures in this article the adiabatic states and their populations are represented by the line type (solid, dashed, etc.) and the colour of the line represents the diabatic states, although of course the latter convention is not entirely rigorous because diabatic states are not unique. The potential energy curves depicted in Fig. 1(b) show “step” at r_NH = 1.45 Å due to the state-averaging over four electronic states, 1¹A₁, 1¹A₂, 1¹B₁, 2¹A₁ (for r_NH < 1.45 Å) and 1¹A₁, 1¹A₂, 1¹B₁, 2¹B₁ (for r_NH > 1.45 Å). When the character of the fourth state changes, the CASSCF optimised orbitals and the CI coefficients change. As a result the calculated energies change slightly as well. However, the fourth electronic excited state is indispensable to reproduce the order of the 3 lower states (1¹A₁, 1¹A₂, 1¹B₁) involved in the dynamics and the dissociative character of ¹πσ* states. The averaging over more than five states is not converged in the right order, at this level of theory (CAS(8,7)-SCF/(aug-)cc-pVDZ). The step is clearly visible under the C_2v symmetry, but it becomes much smoother if the symmetry is lowered. The actual trajectories never reach close proximity of C_2v geometry and they are never affected by “steps”. For this study, the interesting dynamics occurs after passing the two first intersections of the 1¹A₂ and 1¹B₁ states with the state 1¹A₁, and after the departure of the H atom, at r_NH > 2.0 Å. The “step” therefore hardly affects these dynamics. Besides, the N–H cleavage is mostly coupled with the out-of-plane motions of the pyrrolyl ring, which are not under the C_2v symmetry. Hence this level of electronic structure calculation is acceptable.

Table 1 Primary electron configurations, in the SA4-CAS(8,7)-SCF level of theory, of the lowest three electronic states of pyrrole. The equilibrium geometry is under C_2v symmetry, but all the calculations have been performed with no symmetry restrictions. The N–H bond length (r_NH) is 0.996 Å

	σ (9a1)	π (1b1)	π (2b1)	π (1a2)	σ* (10a1)	π* (3b1)	π* (2a2)
1^1A1 (S0)	2	2	2	2	0	0	0
1^1A2 (S1)	2	2	2	1	1	0	0
1^1B1 (S2)	2	2	1	2	1	0	0

---

Table 2 Primary electron configurations, in the SA4-CAS(8,7)-SCF level of theory, of the lowest three electronic states of pyrrole at r_NH = 5.0 Å. The system is no longer molecular pyrrole but a pyrrolyl radical and a hydrogen atom. The lowest two states correspond to the ¹πσ* states while the second excited (1¹A₁) state is correlated with the ¹σσ* state of pyrrole. The geometry depicted below has C_2v symmetry, but electronic structure calculations have been performed with no symmetry restrictions

	π (1b1)	σ (9a1)	π (2b1)	π (1a2)	σ* (10a1)	π* (3b1)	π* (2a2)
1^1A2 (S0)	2	2	2	1	1	0	0
1^1B1 (S1)	2	2	1	2	1	0	0
1^1A1 (S2)	2	1	2	2	1	0	0

---

Fig. 1 (a) One-dimensional potential energy curves of the low-lying electronic states as a function of the NH stretching coordinate. In this figure only the NH bond length was stretched and the molecular frame was frozen. The curves calculated at the SA5-CAS(8,7)-SCF/(aug-)cc-pVDZ level are plotted as solid lines with an offset of 0.8 eV. The 1¹A₁, 1¹A₂, 1¹B₁, and 2¹A₁ states are shown in blue, red, green, and orange, respectively. The dotted lines are from ref. 6. The line colours of the dotted lines correspond to those of the solid lines. The actual dynamics occurs far away from the geometry shown in this figure. (b) Low-lying electronic states 1¹A₁, 1¹A₂, 1¹B₁, 2¹A₁ and 2¹B₁ of pyrrole under C_2v symmetry shown in blue, red, green, orange and black, respectively. The adiabatic states S₀, S₁, S₂, and S₃ are represented by dashed, solid, dotted, and dash-dotted lines.

In the region of r_NH > 2.0 Å, the 1¹A₁ state is no longer the electronic ground state. It finally correlates with the S₂ state, making 1¹A₂ and 1¹B₁ states the ground (S₀) and the first excited (S₁) states respectively. The S₀ and S₁ states correspond to the πσ* states and the S₂ state is the σσ* state. All these surfaces lead to the dissociation of a pyrrole molecule into a pyrrolyl radical and an H atom.

A typical trajectory is depicted in Fig. 2. Fig. 2(a) shows the time evolution of the adiabatic (solid, dashed and dotted lines) and diabatic (labelled with colours) potential energies. Fig. 2(b) depicts the amplitude |d|², which represents the populations of electronic states. The symmetry of the system is not restricted under C_2v so that the electronic states should be labelled as S₀, S₁, S₂; however, these states are strongly correlated with the 1¹A₁, 1¹A₂, 1¹B₁ states of the equilibrium geometry which has C_2v symmetry. In different regions, S₀, S₁, S₂ correspond to different character states 1¹A₁, 1¹A₂, 1¹B₁, which is reflected by the changes in the line colours. At 12.7 fs, the trajectory reaches the S₁–S₀ conical intersection, and almost all the population is transferred from the ¹A₂-like state into the ¹A₁-like state, remaining on the S₁ adiabatic surface. This is in apparent contrast to the behaviour determined from the HRAPTS experiment, where the S₀ (²A₂) state of the radical appears to be populated exclusively. The molecular geometry at the intersection point is determined by a combination of the ν₂₄(b₁) and ν₁₆(b₂) normal modes. The trajectory goes through the second conical intersection between the S₂–S₁ (¹B₁–¹A₁) states at 17.0 fs. Then, roughly 80% of the population of the ¹A₁-like state rapidly moves into the ¹B₁-like state, again remaining on the S₁ surface. The geometry at this second intersection point corresponds to a combination of the ν₁₀(a₂) and ν₁₆(b₂) normal modes. The N–H bond length, shown in Fig. 2(c), is now 1.939 Å, which is twice the equilibrium bond length, so that the pyrrole molecule has essentially already dissociated into the pyrrolyl radical and the hydrogen atom. After passing this conical intersection, the population of S₂, which is now in the ¹A₁-like state, becomes almost constant (|d|² ∼ 0.23). On the other hand, the lowest two S₀ and S₁ states (1¹A₂-like and 1¹B₁-like states of pyrrole) are still coupled. At ∼55 fs, there is a drastic population transfer between the ¹A₂-like and ¹B₁-like state (finally populating the S₀ adiabat) due to another conical intersection, at which the ν₉(a₂), ν₂₀(b₁) and some a₁-symmetry (ν₃, ν₄, ν₇, ν₈, etc.) vibrational modes of the pyrrolyl radical are involved. There are several more significant population exchanges between the S₁–S₀ states of the system (the D₁–D₀ states of the pyrrolyl radical) when the pyrrolyl radical reaches the vicinity of the conical intersection again at 125 and 136 fs.

Fig. 2 (a) Evolution of the potential energies of the S₀, S₁, and S₂ states (solid, dashed and dotted lines) as a function of time during a typical trajectory. The line colours reflect the correlation of adiabatic states with states 1¹A₁, 1¹A₂, and 1¹B₁. These correlations change at the moment of passing in the vicinity of the various intersections, indicated by colour changes. (b) The Ehrenfest populations |d|² of the electronic states S₀, S₁, and S₂ (dashed, solid and dotted lines) along a typical Ehrenfest trajectory. The changes in the line colour reflect the changing correlations with the states 1¹A₁, 1¹A₂, and 1¹B₁, as the trajectory passes through various conical intersections. (c) N–H bond distance as a function of time and its derivative (effectively equivalent to the velocity of dissociating hydrogen atom) along the trajectory shown.

All dissociative Ehrenfest trajectories show similar behaviour. As can be seen from Fig. 2(c) the N–H bond breaks very quickly (within a few tens of femtoseconds) in the vicinity of the intersections around r_(N–H) = 2 Å, shown in Fig. 1. During the H atom departure, the pyrrolyl radical remains in a mixture of two electronic states (with the major component of the population in the S₁ state), which are populated by the dissociative wavepacket as it crosses the initial conical intersections in the r_(N–H) coordinate. However, the average population of the S₁ state (¹B₁-like state of pyrrole) decreases, which is reflected by the fact that any given trajectory spends very little time in this state. Most of the time, the Ehrenfest amplitude of the S₁ state is close to zero and if this state is populated, the trajectory quickly returns to the lower state S₀ (¹A₂-like state of pyrrole). The energy present in the pyrrolyl system is sufficient to allow population of the upper electronic (²B₁) state, yet population of this state is dynamically disfavoured. There appears to be strong coupling between this state and the ground state, such that a statistical population distribution is eventually established and the excess energy tends to be present as nuclear motion (kinetic energy of the wavepacket) on the ground state.

A more traditional way to represent dynamics would simply be to plot the time dependent population of the adiabatic states averaged over all Ehrenfest trajectories, which is shown in Fig. 3. The figure clearly demonstrates two processes. First, the fission of the N–H bond occurs on the time scale of less than 20 fs. Then, the two states of pyrrolyl exchange amplitude and the population of the excited (²B₁) state is transferred to the ground state over a much longer time scale of hundreds of fs. On this time scale the energy is also exchanged between pyrrolyl and hydrogen atom translation motion. Frames 3(a) and (b) show the time dependent population of the electronic states averaged over 3 propagating “bits” obtained for two different choices of the parameter γ in eqn (3), which is the width of the Gaussian Coherent State. For the frame (a), γ = 1 a.u. for all atoms. For the frame (b), the choice suggested in ref. 50 has been employed. Fully converged quantum calculations would be independent of γ. Fig. 3 indicates that the result depends on the choice of the width parameter only slightly.

Fig. 3 Time evolution of populations of the states with different Gaussian Coherent State basis sets. Broken, solid, and dotted lines represent the S₀, S₁, and S₂ states, respectively. Frame (a) is for γ_H/m_H = 5.443 × 10⁻⁴ bohr⁻² emu⁻¹, γ_C/m_C = 4.570 × 10⁻⁵ bohr⁻² emu⁻¹, γ_N/m_N = 3.919 × 10⁻⁵ bohr⁻² emu⁻¹, which corresponds to γ = 1 a.u. for all atoms. Frame (b) is for γ_H/m_H = 5.062 × 10⁻³ bohr⁻² emu⁻¹, γ_C/m_C = 2.076 × 10⁻³ bohr⁻² emu⁻¹, γ_N/m_N = 1.488 × 10⁻³ bohr⁻² emu⁻¹, which was suggested in ref. 50.

The electronic structure of isolated pyrrolyl can give an idea about the electronic structure of pyrrole with strongly elongated N–H bonds. The intersection between the two electronic states of pyrrolyl has been identified previously in ref. 39 and our dynamics calculations also imply its existence. To confirm this, a conical intersection search has been performed at the same level of theory used for the electronic structure in the dynamics calculations. By optimizing both the D₀ (²A₂) minimum and the intersection point, an intrinsic coordinate could be constructed by linear interpolation of the multidimensional coordinate space. By calculation of the state-averaged electronic potentials along this coordinate, the form of the potential and the energetic location of the intersection were established and are displayed in Fig. 4. The intersection is found to occur ∼0.45 eV above the D₀ minimum and is essentially isoenergetic with the D₁ (²B₁) minimum. This D₁ state is therefore coupled at its minimum to the lower state. Due to the presence of this ‘photochemical funnel’, the final population of the electronically excited state is low and the assignment of the strong TKER peaks at around this energy to the high vibrational states of the electronic ground state (X²A₂) rather than to the low vibrational states of the higher electronic state ²B₁ of pyrrolyl is placed on a more solid theoretical footing.

Fig. 4 (a) Cuts through the PESs of the pyrrolyl radical, showing the lowest three states at the SA3-CAS(7,6)-SCF level with (aug)-cc-pVDZ basis. The unit of the intrinsic coordinate is defined as linear displacement of the D₀ minimum from the D₀/D₁ CI geometries (see text). The ground state has a minimum of ²A₂ character, but the D₁ state is intersected at its minimum and there is thus no minimum with ²B₁ character. The states ²A₂, ²B₂ and ²A₁ of pyrrolyl correlate with ¹A₂, ¹B₂ and ¹A₁ of pyrrole with strongly elongated N–H bonds, as the removal of the H atom does not alter the C_2v symmetry. (b) D₀ minimum and (c) D₀/D₁ intersection geometries are shown for reference.

Fig. 5 demonstrates the time dependence of the quantum coupling coefficients D in the wave function (1) for one of the three bundles of 10 Ehrenfest trajectories run in the current simulation. Again frames (a) and (b) represent the result for two different choices of the CS width parameter γ, chosen as γ = 1 a.u. for all atoms (a) and according to the recipe of ref. 50. In both cases quantum coupling is taken into account up to approximately 50 fs. After that, the amplitudes of Ehrenfest trajectories become constant, which indicates that the trajectories have diverged. The CSs, which dress the trajectories, no longer overlap and are therefore no longer coupled. Thus, on the time scale of more than ∼50 fs, the current calculation is only semiclassical, although this has been shown previously to be a good approximation. The time scale of the fully quantum part of the simulation increases with the increase in the size of the trajectory basis. Many more configurations and their trajectories are required for a fully quantum MCE simulation.

Fig. 5 Time dependence plots of the coefficient D. Frame (a) is for γ_H/m_H = 5.443 × 10⁻⁴ bohr⁻² emu⁻¹, γ_C/m_C = 4.570 × 10⁻⁵ bohr⁻² emu⁻¹, γ_N/m_N = 3.919 × 10⁻⁵ bohr⁻² emu⁻¹, which corresponds to γ = 1 a.u. for all atoms. Frame (b) is for γ_H/m_H = 5.062 × 10⁻³ bohr⁻² emu⁻¹, γ_C/m_C = 2.076 × 10⁻³ bohr⁻² emu⁻¹, γ_N/m_N = 1.488 × 10⁻³ bohr⁻² emu⁻¹, which was suggested in ref. 50. Quantum coupling is initially strong and lessens as the trajectories separate. Beyond 50 fs, the amplitudes are constant and the simulation becomes effectively semi-classical.

“On-the-fly” dynamics is computationally expensive. In this particular case we use the SA4-CAS(8,7)-SCF level of electronic structure calculation with the partially augmented cc-pVDZ basis, which requires ∼1600 seconds of CPU time for one time step (in this study the time step was set to 0.1 fs) propagation on an Intel(R) Xeon(R) 5130 processor (2.0 GHz). The advantage of our ab initio MCE approach in its current implementation is that it is a “post-processing” technique. More Ehrenfest trajectories can be run independently and stored. Eqn (7) for their amplitudes D can be run separately after a sufficient number of Ehrenfest trajectories are accumulated. This attractive feature follows from the linear approximation (9) and (10) of the matrix elements. The linear approximation is justified by the fact that in many dimensions only the CSs which are very close to each other overlap and are coupled. In future, the approximation (9) and (10) can be improved by involving more information about the potential energy and its gradient in the neighbouring points of the trajectory. We are currently running more Ehrenfest trajectories and accumulating more data in order to extend the temporal range of the fully quantum part of the simulation and to improve the statistics, which perhaps will allow us to reproduce the translational energy distribution observed in ref. 51. However, even the limited statistics of the current simulation reveals important information about the dissociation dynamics of pyrrole, due to the full dimensionality approach used to simulate the nuclear dynamics.

Thus, we are able to draw strong conclusions about the role of the various conical intersections in the dissociation dynamics of the pyrrole and the subsequent internal dynamics of the nascent pyrrolyl radicals and provide an interpretation for the observation that only the ground (²A₂) state of pyrrolyl appears to be populated at long times in the HRAPTS experiments.

The HRAPTS data indicate exclusive population of the ²A₂ state of pyrrolyl, following N–H bond fission, suggesting that the coupling at the ¹A₂/¹A₁ and ¹A₁/¹B₁ intersections which lie along the dissociation coordinate is minimal and that the dissociative wavepacket remains largely on the diabatic ¹A₂ state. Previous (limited dimensionality) dynamics simulations have predicted population of both the ²A₂ and ²B₁ states, in conflict with observation. Our full dimensionality simulation indicates significant coupling at both intersections, yielding a population of ∼0.8 in the ¹B₁ state after 25 fs (Fig. 3), in agreement with previous dynamics calculations. However, in contrast to the limited dimensionality models, our simulation shows further strong coupling between the D₁ (²B₁) and D₀ (²A₂) states of nearly isolated pyrrolyl, during the dissociation. This leads to decay of the ²B₁ population on the timescale of a few hundred fs and predicts exclusive ²A₂ population on the ns timescale, as observed experimentally. A possible role of this intersection has been proposed in ref. 5 but a multidimensional simulation allows one to demonstrate its importance very naturally.

It is perhaps interesting to consider the timescale of the experiment more carefully. The nascent H atoms show a kinetic energy distribution which is determined by the final pyrrolyl product state distribution. That distribution is only established after a few hundred fs, when the H atom is rather distant (5–7 Å) from the pyrrolyl fragment. Despite the use of nanosecond laser pulses to dissociate the pyrrole, the H atoms depart within 20–30 fs of photon absorption. At this time, the pyrrolyl population is largely on the S₁ (²B₂) surface. One might therefore expect the HRAPTS data to represent a measure of the state of the system after this short time evolution, where after the departure of the H atom the KE distribution is unaffected by the internal dynamics of the pyrrolyl. This appears not to be the case, suggesting that the wavefunctions of the H atom and the partner radical remain coupled even at large spatial separations, as might be expected if H atom is strongly delocalized by preparation of the wave packet and its dynamics.

5. Conclusions

The availability of “on-the-fly” Ehrenfest quantum dynamics allows a more complete analysis of the dissociation of pyrrole than has previously been possible. By including all degrees of freedom in our numerical experiment we do not bias the simulation by the choice of dimensions to a model, hence allowing any and all nuclear coordinates to influence the dynamics. Although the current simulation has been performed mostly semiclassically and has not yet been converged on a fully quantum level, it provides clear evidence that the eventual population of the ²B₁ state of pyrrole is low and thus lends strength to the assignments in the original HRAPTS experiments. It appears that taking into account a large number of degrees of freedom naturally captures the dynamics occurring as a result of several conical intersections, which is harder to see by low dimensionality models or imply from the experimental data alone. It shows that after N–H bond fission the internal dynamics of the pyrrolyl radical ensures coupling between its two lowest electronic states and affects their branching ratio by leaking the population from the upper state to the lower one on the timescale of a few hundred fs. Counterintuitively, the HRAPTS experiment seems sensitive to the final product state distribution, rather than the energy partitioning at the moment of classical dissociation.

Acknowledgements

KS and DS acknowledge the support of EPSRC grant EP/I014500/1.

References

1. J. Wei, J. Riedel, A. Kuczmann, F. Renth and F. Temps, Faraday Discuss., 2004, 127, 267–282.
2. C. T. Middleton, K. de La Harpe, C. Su, Y. K. Law, C. E. Crespo-Hernandez and B. Kohler, Annu. Rev. Phys. Chem., 2009, 60, 217–239.
3. A. L. Sobolewski and W. Domcke, Chem. Phys., 2000, 259, 181–191.
4. A. L. Sobolewski, W. Domcke, C. Dedonder-Lardeux and C. Jouvet, Phys. Chem. Chem. Phys., 2002, 4, 1093–1100.
5. Z. Lan, A. Dupays, V. Vallet, S. Mahapatra and W. Domcke, J. Photochem. Photobiol., A, 2007, 190, 177–189.
6. V. Vallet, Z. G. Lan, S. Mahapatra, A. L. Sobolewski and W. Domcke, J. Chem. Phys., 2005, 123, 144307.
7. M. Barbatti, J. Pittner, M. Pederzoli, U. Werner, R. Mitric, V. Bonacic-Koutecky and H. Lischka, Chem. Phys., 2010, 375, 26–34.
8. V. Vallet, Z. Lan, S. Mahapatra, A. L. Sobolewski and W. Domcke, Faraday Discuss., 2004, 127, 283–293.
9. Z. Lan, W. Domcke, V. Vallet, A. L. Sobolewski and S. Mahapatra, J. Chem. Phys., 2005, 122, 224315.
10. S. Mahapatra, Acc. Chem. Res., 2009, 42, 1004–1015.
11. M. N. R. Ashfold, B. Cronin, A. L. Devine, R. N. Dixon and M. G. D. Nix, Science, 2006, 312, 1637–1640.
12. M. N. R. Ashfold, A. L. Devine, R. N. Dixon, G. A. King, M. G. D. Nix and T. A. A. Oliver, Proc. Natl. Acad. Sci. U. S. A., 2008, 105, 12701–12706.
13. M. N. R. Ashfold, G. A. King, D. Murdock, M. G. D. Nix, T. A. A. Oliver and A. G. Sage, Phys. Chem. Chem. Phys., 2010, 12, 1218–1238.
14. A. L. Devine, B. Cronin, M. G. D. Nix and M. N. R. Ashfold, J. Chem. Phys., 2006, 125, 184302.
15. A. L. Devine, M. G. D. Nix, R. N. Dixon and M. N. R. Ashfold, J. Phys. Chem. A, 2008, 112, 9563–9574.
16. G. A. King, A. L. Devine, M. G. D. Nix, D. E. Kelly and M. N. R. Ashfold, Phys. Chem. Chem. Phys., 2008, 10, 6417–6429.
17. G. A. King, T. A. A. Oliver, M. G. D. Nix and M. N. R. Ashfold, J. Chem. Phys., 2010, 132, 064305.
18. M. G. D. Nix, A. L. Devine, B. Cronin and M. N. R. Ashfold, J. Chem. Phys., 2007, 126.
19. M. G. D. Nix, A. L. Devine, B. Cronin and M. N. R. Ashfold, Phys. Chem. Chem. Phys., 2006, 8, 2610–2618.
20. M. G. D. Nix, A. L. Devine, B. Cronin, R. N. Dixon and M. N. R. Ashfold, J. Chem. Phys., 2006, 125.
21. M. G. D. Nix, A. L. Devine, R. N. Dixon and M. N. R. Ashfold, Chem. Phys. Lett., 2008, 463, 305–308.
22. A. G. Sage, M. G. D. Nix and M. N. R. Ashfold, Chem. Phys., 2008, 347, 300–308.
23. A. J. van den Brom, M. Kapelios, T. N. Kitsopoulos, N. H. Nahler, B. Cronin and M. N. R. Ashfold, Phys. Chem. Chem. Phys., 2005, 7, 892–899.
24. B. Cronin, M. G. D. Nix, R. H. Qadiri and M. N. R. Ashfold, Phys. Chem. Chem. Phys., 2004, 6, 5031–5041.
25. B. Cronin, A. L. Devine, M. G. D. Nix and M. N. R. Ashfold, Phys. Chem. Chem. Phys., 2006, 8, 3440–3445.
26. D. A. Blank, S. W. North and Y. T. Lee, Chem. Phys., 1994, 187, 35–47.
27. H. Yu, N. L. Evans, V. G. Stavros and S. Ullrich, Phys. Chem. Chem. Phys., 2012, 14, 6266–6272.
28. R. Montero, A. P. Conde, V. Ovejas, M. Fernandez-Fernandez, F. Castano, J. R. V. d. Aldana and A. Longarte, J. Chem. Phys., 2012, 137, 064317.
29. K. Saita and D. V. Shalashilin, J. Chem. Phys., 2012, 137, 22A506.
30. D. V. Shalashilin, Faraday Discuss., 2011, 153, 105–116.
31. D. V. Shalashilin, J. Chem. Phys., 2010, 132, 244111.
32. D. V. Shalashilin, J. Chem. Phys., 2009, 130, 244101.
33. B. G. Levine and T. J. Martinez, J. Phys. Chem. A, 2009, 113, 12815–12824.
34. T. J. Martinez, M. Ben-Nun and R. D. Levine, J. Phys. Chem., 1996, 100, 7884–7895.
35. B. G. Levine, J. D. Coe, A. M. Virshup and T. J. Martinez, Chem. Phys., 2008, 347, 3–16.
36. M. Ben-Nun and T. J. Martinez, Chem. Phys., 2000, 259, 237–248.
37. D. Mendive-Tapia, B. Lasorne, G. A. Worth, M. J. Bearpark and M. A. Robb, Phys. Chem. Chem. Phys., 2010, 12, 15725–15733.
38. G. A. Worth, M. A. Robb and B. Lasorne, Mol. Phys., 2008, 106, 2077–2091.
39. A. J. Gianola, T. Ichino, R. L. Hoenigman, S. Kato, V. M. Bierbaum and W. C. Lineberger, J. Phys. Chem. A, 2004, 108, 10326–10335.
40. N. L. Doltsinis and D. Marx, J. Theor. Comput. Chem., 2002, 01, 319–349.
41. P. J. K. H.-J. Werner, G. Knizia, F. R. Manby, M. Schütz, et al., http://www.molpro.net, 2010.
42. R. C. Brown and E. J. Heller, J. Chem. Phys., 1981, 75, 186–188.
43. A. Kirrander and D. V. Shalashilin, Phys. Rev. A, 2011, 84, 13.
44. T. H. Dunning Jr, J. Chem. Phys., 1989, 90, 1007–1023.
45. H. Nakano, T. Tsuneda, T. Hashimoto and K. Hirao, J. Chem. Phys., 1996, 104, 2312–2320.
46. O. Christiansen, J. Gauss, J. F. Stanton and P. Jorgensen, J. Chem. Phys., 1999, 111, 525–537.
47. B. O. Roos, P. A. Malmqvist, V. Molina, L. Serrano-Andres and M. Merchan, J. Chem. Phys., 2002, 116, 7526–7536.
48. P. Celani and H. J. Werner, J. Chem. Phys., 2003, 119, 5044–5057.
49. M. Pastore, C. Angeli and R. Cimiraglia, Chem. Phys. Lett., 2006, 422, 522–528.
50. A. L. Thompson, C. Punwong and T. J. Martínez, Chem. Phys, 2010, 370, 70–77.
51. G. M. Roberts, C. A. Williams, H. Yu, A. S. Chatterley, J. D. Young, S. Ullrich and V. G. Stavros, Faraday Discuss. 10.1039/c2fd20140b.

---