Unveiling anharmonic coupling by means of excited state ab initio dynamics: application to diarylethene photoreactivity

Maria Gabriella Chiariello , Umberto Raucci , Federico Coppola and Nadia Rega *
Dipartimento di Scienze Chimiche, Universita di Napoli Federico II, Napoli, Italy. E-mail: nadia.rega@unina.it

Received 24th July 2018 , Accepted 1st October 2018

First published on 1st October 2018


In this work, excited state ab initio molecular dynamics together with a time resolved vibrational analysis is employed to shed light on the vibrational photoinduced dynamics of a well-known diarylethene molecule experiencing a ring opening reaction upon electronic excitation. The photoreactivity of diarylethenes is recognized to be controlled by a non-adiabatic intersection point between the ground and the first excited state surfaces. The computation of an energy scan, along a suitable reaction coordinate, allows us to identify the region of potential energy surfaces in which the ground (S0) and the first excited (S1) state are well separated. The adiabatic sampling of that region in S1 shows that in the first 3 picoseconds, the central CC bond, which is subject to break, oscillates in an antiphase with respect to the energy gap ΔE(S1 − S0). A multiresolution analysis based on the wavelet transform was then applied to the structural parameters extracted from the excited state dynamics. The wavelet maps show characteristic oscillations of the frequencies, mainly CC stretching and CCC bending localized on the central 4-ring moiety. Moreover, we have identified the main frequency (methyl wagging motion) involved in the modulation of these oscillations. The anharmonic coupling within a group of vibrational modes was therefore highlighted, in good agreement with experimental evidence. For the first time, a quantitative analysis of time resolved signals from a wavelet transform/ab initio molecular dynamics approach was performed.


1 Introduction

Diarylethene molecules are an important class of compounds known to have photochromic properties.1–4 Photochromism is defined as the reversible interconversion, triggered by irradiation, of a chromophore between two or more isomers with different absorption spectra.5,6 Phytochrome and rhodopsin are renowned examples of biological systems showing photochromic activity.7 Aiming to mimic the efficiency of natural biological macromolecules, a great number of synthetic photochromic molecules have been reported so far. These include the azobenzene, spiropyran and diarylethene classes.8–13

The change of the spectroscopic properties is always ascribed to the alteration of the electronic–nuclear structure upon photoexcitation. Therefore, molecules with photochromic activity find many applications as light-driven molecular machines, photoswitching devices or molecular sensors.14–17 Many efforts have been devoted to synthesizing photochromic systems with desired properties, namely: thermal stability of the isomers, resistant properties to side-reactions, sensitivity, and rapid response. In this context, the diarylethene family, due to their excellent performance, has attracted growing attention, more than other classes of photochromic chromophores.18 In principle, photochromic reactions should be fully reversible, and for some diarylethene derivatives, the coloration/decoloration cycles can be repeated more than 105 times.5,18Fig. 1 shows a well-known diarylethene, 1,2-bis(2,4-dimethyl-5-phenyl-3-thienyl)perfluoro-cyclopentene (hereafter named CHD), one of the most used photochromic molecules in practical applications. Upon photoexcitation, CHD undergoes a ring opening reaction, involving the central cyclohexadiene ring, in such a way that photoisomerization converts between the closed and open ring isomer (see Fig. 1). The π conjugation of the open ring isomer is localized on each thiophene ring, while in the closed ring form, the π conjugation is delocalized throughout the molecule, leading to the shift of the absorption spectra to a longer wavelength.18


image file: c8cp04707c-f1.tif
Fig. 1 Open and closed structures of 1,2-bis(2,4-dimethyl-5-phenyl-3-thienyl)perfluoro-cyclopentene (CHD).

Indeed, the closed isomer has an absorption spectrum in the visible (VIS) region (≈600 nm), while the absorption of the open ring species is localized in the ultraviolet (UV) region (≈300 nm). Therefore, photoexcitation with VIS wavelength radiation induces the ring opening reaction, while the opposite ring closing process takes place when the open ring isomer is excited with UV light.18 The photochromic activity of diarylethenes is recognized to be controlled by non-adiabatic decay by means of a conical intersection (CI) seam.19,20 More closely, the ground (S0) and first (S1) excited electronic states become degenerate somewhere along the reaction coordinate and the photoexcited molecular system decays to the ground state in a radiationless manner. CHD was the subject of several experimental studies, where its photochromic properties were exploited for fluorescence signaling and encoding applications.14,15 Moreover, CHD was also deeply investigated in several spectroscopic studies aimed toward elucidating the electronic nuclear dynamics at the origin of its photoreactivity.21–23

In particular, time resolved vibrational spectroscopy (TRVS) experiments revealed characteristic oscillations of some time resolved frequencies.21 Oscillation in the vibrational bands can be ascribed to anharmonic coupling between the high and low frequency modes in the excited state. Quantitative measurement of the anharmonic coupling remains a challenging but fascinating task, because it provides unique information about the shape of the potential energy surfaces (PESs) and the energy exchange between vibrational normal modes.

TRVS experiments corroborate the hypothesis that vibrational modes (mainly stretching and bending motions) localized on the central 4-ring moiety are anharmonically coupled to low frequency methyl wagging modes. These modes are defined as ‘reactive’ because that coupling would lead to the open ring product. Vibrational modes delocalized on the peripheral phenyl groups are instead labeled as unreactive, leading the system back to the ground state.

Of course, unveiling at the molecular level the driving forces controlling the CHD photoreactivity is extremely important to plan the rational design of a photochromic system with better performance.

For this reason, in this paper, we aim to study the electronic–nuclear driving forces that lead the photoexcited system toward the conical intersection region. In particular, we focus on the analysis of the vibrational modes responsible for the photoreactive process, and the coupling between them. Theoretical studies of CHD are quite unusual, due to the large size of the system. In principle, non-adiabatic photoreactivity should be described by Multi-Configurational Self Consitent Field (MC-SCF) methods, allowing for a suitable characterization of the PES region where the crossing between the two electronic surfaces takes place.24–28 Unfortunately, they are not feasible for large size systems such as CHD due to their high computational cost.

On the other hand, the regions of the PES(s) in which the ground and excited state surfaces are well separated (namely regions in which the Born–Oppenheimer approximation still holds) can be accurately investigated by methods rooted in Density Functional Theory (DFT) and its time dependent version (TD-DFT).8,29–42

Indeed, these regions far from the conical intersection branching space can be accurately described within the adiabatic approximation by running single surface excited state dynamics and getting information about the evolution of the vibrational properties by means of a rational analysis of the trajectories.

The use of ab initio molecular dynamics (AIMD) methods is here associated with a time resolved vibrational analysis based on the wavelet transform.43–45 The application of the wavelet transform combined with excited state AIMD was already tested and validated, giving promising results in studies aimed at disentangling the time evolution of a simulated Stokes shift,46 and of the dipole moment in exciton dynamics.47 AIMD combined with wavelet analysis was also successfully adopted to describe the photoreactivity of green fluorescent protein in the time–frequency domain, and to simulate spectroscopic signals of the pyranine photoacid in water solution.36,37

Following this line of research, the vibrational analysis is here expanded to investigate anharmonic coupling between vibrational signals of the CHD photoswitch. More closely, the time resolved vibrational bands, as a result of the wavelet transform, are further analyzed with the purpose of recognizing the frequency ruling the characteristic oscillations of time resolved vibrational signals. The wavelet maps show characteristic oscillations of the frequencies, mainly CC stretching and CCC bending localized on the central 4-ring moiety. Moreover, we have identified the main frequency (methyl wagging) involved in the modulation of these oscillations. The anharmonic coupling within this group of vibrational modes was, therefore, highlighted, in good agreement with experimental evidence.

This paper is organized as follows: after a brief description of the computational details and the wavelet based vibrational analysis (Section 2), the static characterization of the reaction profiles is discussed in Section 3.1. The dynamical sampling and the vibrational analysis are then presented in Sections 3.2 and 3.3. Final remarks and future perspectives are given as conclusions.

2 Models and methods

2.1 Simulation details

In order to characterize the ground and the excited state PESs of the CHD closed ring isomer (see Fig. 1), DFT and TD-DFT levels of theory were employed. The global hybrid functional B3LYP was adopted in combination with the 6-31g(d,p) basis set.48–50 This method allowed us to reproduce the spectroscopic features of the chromophore, obtaining photophysical signatures in fair agreement with the experimental absorption and fluorescence spectra (see Section 3.1 for more details).

Starting from the CHD minimum energy structure, constrained energy profiles were built by partial relaxed structural optimizations along the central CC distance, hereafter CC1, representing a reasonable simplified reaction coordinate for the ring opening reaction (see Fig. 1). The CC1 distance was scanned from 1.54 Å to 3.98 Å, which correspond to the values that the CC1 bond assumes, respectively, in the closed and open ring isomers in their ground state minimum energy structures. The intermediate points were obtained by keeping the CC1 distance fixed while relaxing in the ground state all the other degrees of freedom. The S1 excited state profile was vertically computed in these points at the TD-B3LYP/6-31g(d,p) level of theory.

Because of the non-adiabatic photoreactivity of the CHD molecule, a quantitative measure of the coupling between the ground and first excited state surfaces is helpful. Therefore, the ground-to-excited state non-adiabatic coupling (NAC) was computed along the points of the energetic profile. In particular, the first order non-adiabatic coupling matrix element can be written as:

 
image file: c8cp04707c-t1.tif(1)
where En and Em are the energy eigenvalues of the states Ψn and Ψm, [V with combining circumflex]ne denotes the operator of the electron–nuclear interaction and ξ indicates a Cartesian nuclear coordinate. In this work, the first order NAC was computed analytically in the framework of linear response TD-DFT,51,52 namely the excited state wavefunction is not computed, since all information about the excited state properties is contained in the time evolution of the ground reference state. Regarding the ab initio molecular dynamics simulation, we used the same level of theory of static characterization, namely B3LYP/6-31g(d,p). The ground state sampling was performed by means of the Atom-Centered Density Matrix Propagation (ADMP) approach,53–61 collecting a trajectory for 5 ps with a time step of 0.2 fs and maintaining a constant temperature of 300 K. Starting from the sampling of the ground state phase space, a proper number of points (positions and momenta) should be chosen in order to represent, on average, the thermal distribution of the nuclear structure and the kinetic energy. Due to the number of degrees of freedom involved, however, this procedure would imply an unaffordable computational cost. We therefore focused on the sampling of the CC1 distance value: we selected three points with a CC1 distance slightly larger or close to the average value in S0, and with a total kinetic energy close to the average value as well. These configurations were used as the starting points for the propagation of the excited state trajectories. Simulations in S1 were performed following a Born–Oppenheimer AIMD (BOMD)/TD-DFT scheme, in which energies and energy derivatives are computed on-the-fly at each step of the dynamics.62,63 Trajectories were propagated for 3 ps using a time step of 0.7 fs. In all the cases, the energy was conserved with fluctuations within 0.0004 Hartree. All the calculations were carried out with the Gaussian16 suite of programs.64

2.2 Time resolved frequency multiresolution analysis

Upon electronic excitation, the CHD molecule undergoes a complex vibrational dynamics. In this work, the composition of vibrational modes is approximated with structural parameters extracted from the trajectories. They contain the key information about the relaxation dynamics that the modes experience in the excited state.36,37,46 The multiresolution vibrational analysis is based on the wavelet transform, which was already and successfully used in combination with AIMD.47,65,66 The theoretical framework of the wavelet transform is here only briefly reported. The continuous wavelet transform is employed:
 
image file: c8cp04707c-t2.tif(2)

In this case, C(t) is the time dependent structural parameter extracted from the excited state trajectories. ψa,b is member of the so-called wavelet basis, with functions that are the translated and contracted/dilated version of a mother wavelet. We choose the Morlet function as the mother wavelet.67–69 The dilatation/contraction of the mother function is ruled by the a parameter, which is proportional to the inverse of a frequency. The b parameter controls, instead, the translation of the wavelet basis along the temporal axis. We plot the magnitude square of the transform |W(v,t)|2 as the intensity of the instantaneous frequency contribution to the signal. The result of the wavelet integration is a signal localized in both the time and frequency domains. The application of the wavelet allows for investigating the characteristic oscillation behavior of the frequencies resolved in time. The analysis was performed on the three excited state trajectories. In the next section, we present the wavelet transformed signals averaged on the three trajectories. The wavelet analysis carried out on single trajectories is instead reported in the ESI.

3 Results and discussion

3.1 TD-DFT picture of the potential energy surfaces

In the following, the focus will be on the direct ring opening reaction, namely the photo-isomerization starting from the closed ring reagent that leads to the open ring product. The vertical excitation energy (VEE) is obtained from the ground state minimum geometry of the CHD closed ring form. The computed VEE value is 2.12 eV in the gas phase, which is in nice agreement with the value obtained from the experimental absorption spectrum in cyclohexane of 2.21 eV. In principle, the computational/experimental discrepancy could be reduced by including the effect of the solvent. Nevertheless, the difference of 0.09 eV, corresponding to a small gap of 2 kcal mol−1, can be considered negligible to our aims, as well as the effect of the cyclohexane solvent.

The transition leading to the S1 excited state has a π–π* nature, involving frontier orbitals shown in Fig. 2.


image file: c8cp04707c-f2.tif
Fig. 2 HOMO and LUMO contour plots computed for CHD in the ground-state minimum energy structure obtained at the B3LYP/6-31g(d,p) level of theory.

Upon excitation, the main part of the charge redistribution concerns the π conjugation localized on the central rings, namely the thiophene rings, the central cyclohexadiene and the fluorinated cyclopentane, the latter acting as an electron acceptor. On the other hand, the peripheral phenyl rings are only slightly involved in the charge rearrangement associated with the excitation. A side perspective of the orbitals (see lower panel of Fig. 2) makes it possible to observe the new electronic arrangement on the central CC1 bond, i.e. the bond subject to break in the excited state. The electronic charge is strongly reduced along the CC1 bond when going from the ground to the first excited state.

In order to get a first picture of the potential energy surfaces involved in the ring opening reaction, constrained potential energy profiles along a suitable reaction coordinate (namely the CC1 central bond) were calculated in the ground state and in the first excited state (see Section 2.1 for details). In principle, DFT and its time dependent version, as a single reference electronic structure theories, are not able to characterize the crossing points between two electronic surfaces. Nevertheless, this procedure aims at understanding where the two surfaces start to become degenerate along the reaction coordinate, so that it will be possible to identify some regions of the PESs that can be treated in an adiabatic way. The ground and the excited state energy profiles are shown in Fig. 3. The profiles represent relative energies with respect to the energy of the closed ring CHD, which was instead set to zero. Moreover, the ground-to-excited state non-adiabatic couplings are also shown in Fig. 3 (green dashed line). The NAC computation is here employed to report the coupling between the ground and the first excited states along the constrained PES, with the CC1 bond as the reaction coordinate.


image file: c8cp04707c-f3.tif
Fig. 3 Ground (blue) and excited (red) constrained energy profiles calculated for the ring opening reaction at the B3LYP/6-31g(d,p) level of theory. Blue circles indicate energy values of the fully relaxed open and closed ring structures. The non-adiabatic couplings (green dashed line) computed along the profile are also shown. The grey surface represents the non-adiabatic region. The violet area defines the adiabatic region considered in this work.

In the ground state, the transformation between the closed and open ring forms would require overcoming an energy barrier of about 50 kcal mol−1. Instead, when the molecule is photoexcited, the scenario dramatically changes. A small barrier of 6 kcal mol−1 separates the Franck–Condon region from the non-adiabatic zone, where the energy levels of the ground and the first excited states are very close (highlighted by the grey area in Fig. 3). The excited system drops down to the ground state through the conical intersection. The radiationless decay to the ground state leads to the open ring isomer, which is about 13 kcal mol−1 more stable than the closed one. This picture is in fair agreement with MC-SCF studies on diarylethene model systems, which consist of the central 4-ring reactive core. The CI point was identified at a CC1 distance of about 2 Å.70

As already mentioned above, DFT does not allow characterization of the intersection points, although we can get information about the region of the PESs in which the ground and excited state surfaces are well separated and the adiabatic approximation still holds. More closely, within the range of values of the CC1 bond from 1.54 to 1.75 Å, the two electronic surfaces are clearly non degenerate (see the violet area in Fig. 3). Within the given level of theory, the two surfaces become too close in energy at a CC1 distance of about 1.80 Å, where the NAC profile reaches its maximum value. After that, it decreases with increasing CC1 distance. On the other hand, within the 1.54–1.75 Å range, NACs undergo little variations and are always low. Interestingly, in the 1.50–1.75 Å range, we can also observe a significant decrease of the profile slope when passing from S0 to S1, suggesting that the CC1 stretching motion becomes less stiff upon the electronic excitation.

The adiabatic dynamical sampling of that region will be the subject of the next section.

3.2 Excited state ab initio molecular dynamics

A first insight into the response of the CHD nuclei to the photoexcitation and to the new electronic arrangement can be obtained by analyzing the average values of some key CC distances of the central moiety, in both ground and excited states. The average values and labels of all the bonds of the central 4-ring system are shown in Fig. 4. It is worth comparing the ground and excited state behaviors by considering the contours of the frontiers orbitals (HOMO–LUMO in Fig. 2) involved in the electronic transition. When going from the ground to the excited state, some CC distances of the cyclohexadiene ring (CC5 and CC6) and of the thiophene ring (CC9 and CC10) become more stretched, going from a typical double bond value in the ground state (1.37 Å) to the average values of 1.42 (CC5, CC6) and 1.40 Å (CC9, CC10) in the excited state. This change occurs according to a depletion of electron density on these bonds, as is possible to observe from the frontier orbital contours. On the other hand, an increase of density on CC2, CC3, CC4 (cyclohexadiene), CC7, and CC8 (thiophenes) bonds is observed, making the distance shorter in the excited state.
image file: c8cp04707c-f4.tif
Fig. 4 Average distances (Å) obtained from the AIMD simulation of CHD in the ground state (central structure) and the S1 excited state (structure on the right). Structure on the left shows numbering of the bonds.

It is interesting to note that the excited state nuclear rearrangement leads the CC bond in the CC7–CC5–CC2–CC6–CC8 chain to assume the same length of 1.42 Å on average, in the middle of a single and double bond value.

The depletion of electron density on the CC1 bond makes it weaker and longer in the excited state, showing an average value of 1.55 Å in S0 and 1.60 Å in S1. Hence, the CC1 bond shows on average longer distances in S1, as a result of the new electronic arrangement. Fig. 5 shows the time evolution of the CC1 distance in both S0 and S1. Larger oscillations are observed in the excited state, with a maximum value of 1.75 Å. As already shown in the previous section, a CC1 distance of 1.75 Å represents a threshold value for which the non-adiabatic coupling remains low. Therefore, starting from the Franck–Condon region, our trajectories describe the pathway toward the CI zone, although remaining in the adiabatic region. Indeed, oscillations of the CC1 bond strongly affect the S0–S1 energy gap.


image file: c8cp04707c-f5.tif
Fig. 5 (a) Ground (blue) and excited (red) state time evolution of the CC1 distance. (b) Comparison between the normalized fluctuactions of the CC1 distance in the excited state (red line) and the S0 − S1 energy gap ΔE (blue line).

We define the normalized fluctuations of a time evolving function f(t) as

 
image file: c8cp04707c-t3.tif(3)
where [f with combining macron](t) is the average value over time. In Fig. 5, we show normalized fluctuations of the CC1 distance in S1 (red line) and of the S1 − S0 energy difference (ΔES1−S0) (blue line).

The graph clearly shows that the two curves oscillate out of phase, i.e. maxima of the energy difference correspond to minima of the CC1 distance over time. This means that the larger the CC1 distance, the smaller the energy gap between S1 and S0. Moreover, analyzing the behaviour of ΔES1−S0 over time, we computed an average value of 1.50 eV (corresponding to 35 kcal mol−1) with a standard deviation of 0.17 eV, with maximum and minimum values of 2.14 and 0.97 eV, respectively.

3.3 Time resolved vibrational analysis

Upon excitation, CHD undergoes a strong vibrational activity, mirroring the nuclear evolution on the excited state PES. The photoinduced vibrational dynamics was recently analyzed by two-dimensional femtosecond stimulated Raman techniques,21 which were specifically employed to reveal the anharmonic coupling between vibrational modes. In particular, in a recent work,71 authors inferred a quantitative relationship between high frequency oscillating bands and low frequency modes, modulating them by means of an anharmonic coupling constant. In the case of the CHD molecule, the dynamics of a group of vibrational bands was identified as strongly coupled. These bands were labeled as ‘reactive’, because the coupling was responsible for the barrier crossing and the subsequent ring opening reaction through the conical intersection. In particular, the coupling was proven by observing the oscillation patterns in both the intensity and the frequency time evolution of high and low frequency modes.

According to this hypothesis based on experimental data, the reactive vibrational bands include CC stretching modes (at 1181 and 1333 cm−1) and a CCC bending mode (at 467 cm−1), localized on the cyclohexadiene and thiophene rings. These are coupled to the low frequency wagging motions of the methyl groups on the cyclohexadiene ring (at 191 cm−1). Therefore, the photoinduced relaxation of the stretching and bending vibrational bands shows an oscillating behavior modulated by the low frequency methyl wagging mode.21

In the following discussion, the focus is on the modes belonging to the ‘reactive’ group and the coupling between them. First, we identified reactive modes from the frequency calculation of the CHD S1 minimum energy structure, by using the B3LYP/6-31g(d,p) level of theory. Then, we chose the key structural parameters mainly involved in these normal modes and constructed the corresponding averaged relaxation functions from the collected excited state trajectories. Finally, the adopted strategy included a time resolved vibrational analysis based on the wavelet transform of these relaxation functions.

In Table 1, we show the harmonic frequencies calculated for the reactive modes in the CHD S1 minimum energy structure. These values are compared with the corresponding frequencies averaged on the three excited state trajectories. The experimental data are also shown.21 Both simulated and experimental data refer to averages over a period of 3 ps, namely approximately the relaxation time spent on the S1 PES, prior to entering the conical intersection region. Moreover, we recall that the simulated frequencies are intrinsically anharmonic. As a matter of fact, we observe a nice agreement between simulated and experimental values.

Table 1 AIMD and harmonic frequencies (cm−1) of the vibrational modes analyzed in the present work. Harmonic frequencies refer to the S1 minimum energy structure computed at the B3LYP/6-31g(d,p) level of theory. Experimental values from ref. 21
Exp. AIMD Harmonic
Methyl wagging 191 208 214
CC3 stretching 1181 1186 1194
CC7 stretching 1333 1380 1407
CC5–CC7 bending 467 469 509


The AIMD simulation allows for an exploration of the ground and excited state potential energy surfaces that are intrinsically anharmonic, hence a quantitative analysis of the time resolved oscillating frequencies allows one to get information about the anharmonic coupling.

The CC1 distance oscillations, discussed in the previous section, affect other nuclear motions, in particular, the vibrational modes involving the methyl groups of the CC1 bond. In the closed ring form, the methyl groups are basically perpendicular to the thiophene ring plane. The open ring product is instead characterized by the methyl groups in the plane of the thiophene rings. The methyl wagging normal mode at 191 cm−1[thin space (1/6-em)]21 (see the ESI) describes the motion from the perpendicular toward the in plane arrangement, and vice versa. The structural parameter chosen as representative of the methyl wagging mode is the (H3)C–CC1–C(H3) dihedral angle, and the 2D wavelet map of the methyl wagging extracted from the excited state trajectories is shown in Fig. 6. The time resolved vibrational band of the wagging mode is centered at 200 cm−1, in fair agreement with the experimental one.


image file: c8cp04707c-f6.tif
Fig. 6 Left panel: CHD structure with the atoms participating in the methyl wag (AIMD frequency 208 cm−1) highlighted in red. Right panel: 2D wavelet map of the methyl wagging motion in the excited state. The color scale represents the intensity in arbitrary units.

The methyl wagging is the low frequency mode recognized as coupled and modulating the higher frequency modes localized on the central rings. These high frequency modes are the CC stretching on the cyclohexadiene ring (experimental frequency of 1181 cm−1),21 the CC stretching localized on the thiophene rings (exp. 1333 cm−1) and the CCC bending of cyclohexadiene and thiophene rings (exp. 467 cm−1). The corresponding wavelet 2D maps are reported in Fig. 7, 8 and 9, respectively. Two kinds of CC stretching are analyzed: the CC stretching on the central ring and on the thiophene rings, with AIMD frequencies of 1186 cm−1 and 1380 cm−1, respectively.


image file: c8cp04707c-f7.tif
Fig. 7 Left panel: CHD structure with the atoms participating in the CC stretching (AIMD frequency 1186 cm−1) highlighted in red. Central panel: 2D wavelet map of the CC stretching localized on the cyclohexadiene moiety. The color scale represents the intensity in arbitrary units. Right panel: Discrete Fourier transform of the intensity oscillations.

image file: c8cp04707c-f8.tif
Fig. 8 Left panel: CHD structure with the atoms participating in the CC stretching (AIMD frequency 1380 cm−1) highlighted in red. Central panel: 2D wavelet map of the CC stretching localized on the thiophene rings. The color scale represents the intensity in arbitrary units. Right panel: Discrete Fourier transform of the intensity oscillations.

image file: c8cp04707c-f9.tif
Fig. 9 Left panel: CHD structure with the atoms participating in the CCC bending (AIMD Frequency 469 cm−1) highlighted in red. Central panel: 2D wavelet map of the CCC bending localized on the cyclohexadiene and thiophene rings. Right panel: Discrete Fourier transform of the intensity oscillations.

The wavelet maps of the two stretching modes (see Fig. 7 and 8) show clear intensity and frequency oscillations over time. In order to unveil the coupling with the methyl wagging mode, a quantitative analysis of the time dependent frequencies was performed. In particular, frequencies corresponding to the maximum intensity were extracted from the time resolved wavelet output, and a discrete Fourier transform (FT) was performed on the maximum intensity value in order to reveal the oscillation frequency. The discrete Fourier transform of the time resolved maximum intensity associated with the CC stretching modes is shown in Fig. 7 and 8. In both cases, the FT shows that the main frequency modulating the intensity oscillations of the CC stretching mode is at 190 cm−1, matching the methyl wagging frequency.

The 2D wavelet map of the CCC bending of the cyclohexadiene and thiophene rings is quite complex (Fig. 9). The signal of interest is centered at a frequency of 469 cm−1 (AIMD frequency), although other contributions can be observed. As a matter of fact, the CCC bending involving cyclohexadiene and thiophene rings is part of another deformation mode at 690 cm−1.

The wavelet map of the CCC bending also shows a contribution at about 200 cm−1. Because the CCC bending does not participate in any modes around 200 cm−1, the signal appears as a result of the anharmonic coupling with the methyl wags. Moreover, the intensity oscillation corresponding to the frequency centered at 469 cm−1 was extracted and analyzed. The FT of the maximum intensity over time (Fig. 9) shows a signal at 191 cm−1, confirming the coupling with the methyl wagging mode.

Our findings are in agreement with the experimental evidence of a coupling between the methyl wagging (acting as modulating mode) and the CC stretching and CCC bending modes. A similar behavior of low/high frequency anharmonic coupling has been recognized for the photorelaxation processes of the pyranine photoacid37,72 and the photodynamics of the green fluorescent protein, with the subsequent excited state proton transfer reaction.36,73 However, a quantitative analysis of the time resolved bands from AIMD has been applied here for the first time and it has provided satisfactory results.

4 Conclusions

In this work, we studied the photoinduced vibrational dynamics of a diarylethene photoswitch, combining ab initio molecular dynamics simulations and a time resolved vibrational analysis. As a response to electronic excitation, the central CC bond experiences wider oscillations in the sampled time. These fluctuations are out of phase with respect to the energy gap S1 − S0, suggesting that the CC bond is one of the main coordinates ruling the conical intersection. A great deal of information on the process was obtained by analyzing the relevant structural parameters in both the time and frequency domains, by means of the wavelet transform. The wavelet maps revealed characteristic oscillations of the CC stretching and CCC bending modes localized on the central 4-ring diarylethene core. A quantitative analysis of the oscillating behavior allowed us to recognize the frequency that modulates the CC stretching and CCC bending oscillations. Our results indicate that the methyl wagging mode is mainly involved in the modulation of the CC stretching and CCC bending frequencies, revealing a strong coupling between these vibrational modes, as already experimentally suggested.

In conclusion, ab initio molecular dynamics combined with time resolved vibrational analysis appears to be a suitable tool to disentangle nuclear photoinduced dynamics, even in the case of challenging and complex systems like the CHD molecule. The experimental evidence is well reproduced, suggesting the accuracy and the general reliability of the method.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors gratefully acknowledge funding from Gaussian, Inc. (Wallingford, CT).

Notes and references

  1. H. Tian and S. Yang, Chem. Soc. Rev., 2004, 33, 85–97 RSC.
  2. M. Irie, T. Fukaminato, K. Matsuda and S. Kobatake, Chem. Rev., 2014, 114, 12174–12277 CrossRef CAS PubMed.
  3. M. Morimoto, T. Sumi and M. Irie, Materials, 2017, 10, 1021 CrossRef PubMed.
  4. T. Ichikawa, M. Morimoto and M. Irie, Dyes Pigm., 2017, 137, 214–220 CrossRef CAS.
  5. M. Irie, Chem. Rev., 2000, 100, 1683–1684 CrossRef CAS PubMed.
  6. M. Irie, K. Sakemura, M. Okinaka and K. Uchida, J. Org. Chem., 1995, 60, 8305–8309 CrossRef CAS.
  7. S. H. Bhoo, T. Hirano, H.-Y. Jeong, J.-G. Lee, M. Furuya and P.-S. Song, J. Am. Chem. Soc., 1997, 119, 11717–11718 CrossRef CAS.
  8. M. Savarese, U. Raucci, P. A. Netti, C. Adamo, N. Rega and I. Ciofini, Theor. Chem. Acc., 2016, 135, 211–217 Search PubMed.
  9. M. Rini, A.-K. Holm, E. T. Nibbering and H. Fidder, J. Am. Chem. Soc., 2003, 125, 3028–3034 CrossRef CAS.
  10. N. P. Ernsting, Chem. Phys. Lett., 1989, 159, 526–531 CrossRef CAS.
  11. A. A. Beharry and G. A. Woolley, Chem. Soc. Rev., 2011, 40, 4422–4437 RSC.
  12. M. Dong, A. Babalhavaeji, S. Samanta, A. A. Beharry and G. A. Woolley, Acc. Chem. Res., 2015, 48, 2662–2670 CrossRef CAS.
  13. K. Matsuda and M. Irie, J. Photochem. Photobiol., C, 2004, 5, 169–182 CrossRef CAS.
  14. C. Yun, J. You, J. Kim, J. Huh and E. Kim, J. Photochem. Photobiol., C, 2009, 10, 111–129 CrossRef CAS.
  15. M. Morimoto and M. Irie, J. Am. Chem. Soc., 2010, 132, 14172–14178 CrossRef CAS.
  16. T. Fukaminato, T. Sasaki, T. Kawai, N. Tamai and M. Irie, J. Am. Chem. Soc., 2004, 126, 14843–14849 CrossRef CAS.
  17. M. Estrader, J. Salinas Uber, L. A. Barrios, J. Garcia, P. Lloyd-Williams, O. Roubeau, S. J. Teat and G. Arom, Angew. Chem., Int. Ed., 2017, 56, 15622–15627 CrossRef CAS.
  18. M. Irie, Photochem. Photobiol. Sci., 2010, 9, 1535–1542 RSC.
  19. M. Boggio-Pasqua, M. Ravaglia, M. J. Bearpark, M. Garavelli and M. A. Robb, J. Phys. Chem. A, 2003, 107, 11139–11152 CrossRef CAS.
  20. C. L. Ward and C. G. Elles, J. Phys. Chem. A, 2014, 118, 10011–10019 CrossRef CAS PubMed.
  21. D. T. Valley, D. P. Hoffman and R. A. Mathies, Phys. Chem. Chem. Phys., 2015, 17, 9231–9240 RSC.
  22. Y. Ishibashi, T. Umesato, M. Fujiwara, K. Une, Y. Yoneda, H. Sotome, T. Katayama, S. Kobatake, T. Asahi, M. Irie and H. Miyasaka, J. Phys. Chem. C, 2016, 120, 1170–1177 CrossRef CAS.
  23. H. Sotome, T. Nagasaka, K. Une, C. Okui, Y. Ishibashi, K. Kamada, S. Kobatake, M. Irie and H. Miyasaka, J. Phys. Chem. Lett., 2017, 8, 3272–3276 CrossRef CAS PubMed.
  24. B. O. Roos, Adv. Chem. Phys., 2007, 69, 399–445 CrossRef.
  25. B. O. Roos, P. R. Taylor and P. E. Siegbahn, Chem. Phys., 1980, 48, 157–173 CrossRef CAS.
  26. M. Ben-Nun and T. J. MartÌAnez, Chem. Phys., 2000, 259, 237–248 CrossRef CAS.
  27. W. Domcke and D. R. Yarkony, Annu. Rev. Phys. Chem., 2012, 63, 325–352 CrossRef CAS.
  28. L. S. Cederbaum, E. Gindensperger and I. Burghardt, Phys. Rev. Lett., 2005, 94, 113003 CrossRef.
  29. R. E. Stratmann, G. E. Scuseria and M. J. Frisch, J. Chem. Phys., 1998, 109, 8218–8224 CrossRef CAS.
  30. M. E. Casida and M. Huix-Rotllant, Annu. Rev. Phys. Chem., 2012, 63, 287–323 CrossRef CAS PubMed.
  31. M. E. Casida and T. A. Wesołowski, Int. J. Quantum Chem., 2004, 96, 577–588 CrossRef CAS.
  32. C. Adamo, T. Le Bahers, M. Savarese, L. Wilbraham, G. Garca, R. Fukuda, M. Ehara, N. Rega and I. Ciofini, Coord. Chem. Rev., 2015, 304, 166–178 CrossRef.
  33. E. Runge and E. K. Gross, Phys. Rev. Lett., 1984, 52, 997 CrossRef CAS.
  34. A. Petrone, D. B. Lingerfelt, D. B. Williams-Young and X. Li, J. Phys. Chem. Lett., 2016, 7, 4501–4508 CrossRef CAS.
  35. U. Raucci, M. Savarese, C. Adamo, I. Ciofini and N. Rega, J. Phys. Chem. B, 2015, 119, 2650–2657 CrossRef CAS.
  36. G. Donati, A. Petrone, P. Caruso and N. Rega, Chem. Sci., 2018, 9, 1126–1135 RSC.
  37. M. G. Chiariello and N. Rega, J. Phys. Chem. A, 2018, 122, 2884–2893 CrossRef CAS.
  38. P. Cimino, U. Raucci, G. Donati, M. G. Chiariello, M. Schiazza, F. Coppola and N. Rega, Theor. Chem. Acc., 2016, 135, 1–12 Search PubMed.
  39. M. Savarese, U. Raucci, R. Fukuda, C. Adamo, M. Ehara, N. Rega and I. Ciofini, J. Comput. Chem., 2017, 38, 1084–1092 CrossRef CAS.
  40. M. Savarese, P. A. Netti, N. Rega, C. Adamo and I. Ciofini, Phys. Chem. Chem. Phys., 2014, 16, 8661–8666 RSC.
  41. F. Perrella, U. Raucci, M. G. Chiariello, M. Chino, O. Maglio, A. Lombardi and N. Rega, Biopolymers, 2018 DOI:10.1002/bip.23225.
  42. A. Petrone, J. Cerezo, F. J. A. Ferrer, G. Donati, R. Improta, N. Rega and F. Santoro, J. Phys. Chem. A, 2015, 119, 5426–5438 CrossRef CAS PubMed.
  43. C. Torrence and G. P. Compo, Bull. Am. Meteorol. Soc., 1998, 79, 61–78 CrossRef.
  44. M. Farge, N. Kevlahan, V. Perrier and K. Schneider, Wavelets in Physics, 1999, pp. 117–200 Search PubMed.
  45. M. Farge, Annu. Rev. Fluid Mech., 1992, 24, 395–458 CrossRef.
  46. A. Petrone, G. Donati, P. Caruso and N. Rega, J. Am. Chem. Soc., 2014, 136, 14866–14874 CrossRef CAS PubMed.
  47. G. Donati, D. B. Lingerfelt, A. Petrone, N. Rega and X. Li, J. Phys. Chem. A, 2016, 120, 7255–7261 CrossRef CAS.
  48. A. D. Becke, J. Chem. Phys., 1992, 96, 2155–2160 CrossRef CAS.
  49. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098 CrossRef CAS.
  50. A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377 CrossRef CAS.
  51. R. Send and F. Furche, J. Chem. Phys., 2010, 132, 044107 CrossRef PubMed.
  52. D. B. Lingerfelt, D. B. Williams-Young, A. Petrone and X. Li, J. Chem. Theory Comput., 2016, 12, 935–945 CrossRef CAS.
  53. H. B. Schlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, A. D. Daniels, G. E. Scuseria and M. J. Frisch, J. Chem. Phys., 2001, 114, 9758–9763 CrossRef CAS.
  54. S. S. Iyengar, H. B. Schlegel, J. M. Millam, G. A. Voth, G. E. Scuseria and M. J. Frisch, J. Chem. Phys., 2001, 115, 10291–10302 CrossRef CAS.
  55. H. B. Schlegel, S. S. Iyengar, X. Li, J. M. Millam, G. A. Voth, G. E. Scuseria and M. J. Frisch, J. Chem. Phys., 2002, 117, 8694–8704 CrossRef CAS.
  56. S. S. Iyengar, H. B. Schlegel, G. A. Voth, J. M. Millam, G. E. Scuseria and M. J. Frisch, Isr. J. Chem., 2002, 42, 191–202 CrossRef CAS.
  57. N. Rega, S. S. Iyengar, G. A. Voth, H. B. Schlegel, T. Vreven and M. J. Frisch, J. Phys. Chem. B, 2004, 108, 4210–4220 CrossRef CAS.
  58. N. Rega, Theor. Chem. Acc., 2006, 116, 347–354 Search PubMed.
  59. G. Brancato, N. Rega and V. Barone, J. Chem. Phys., 2008, 128, 04B607 CrossRef.
  60. N. Rega, G. Brancato and V. Barone, Chem. Phys. Lett., 2006, 422, 367–371 CrossRef CAS.
  61. D. Branduardi, M. De Vivo, N. Rega, V. Barone and A. Cavalli, J. Chem. Theory Comput., 2011, 7, 539–543 CrossRef CAS.
  62. R. N. Barnett and U. Landman, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 2081 CrossRef CAS.
  63. I. Tavernelli, U. F. Röhrig and U. Rothlisberger, Mol. Phys., 2005, 103, 963–981 CrossRef CAS.
  64. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, et al., Gaussian 16 Revision A.03, Gaussian Inc., Wallingford CT, 2016 Search PubMed.
  65. M. Pagliai, F. Muniz-Miranda, G. Cardini, R. Righini and V. Schettino, J. Mol. Struct., 2011, 993, 438–442 CrossRef CAS.
  66. F. Muniz-Miranda, M. Pagliai, G. Cardini and V. Schettino, J. Chem. Theory Comput., 2011, 7, 1109–1118 CrossRef CAS.
  67. I. Daubechies, IEEE Trans. Inf. Theory, 1990, 36, 961–1005 CrossRef.
  68. O. Rioul and M. Vetterli, IEEE Signal Process. Mag., 1991, 8, 14–38 Search PubMed.
  69. H. Weng and K. Lau, J. Atmos. Sci., 1994, 51, 2523–2541 CrossRef.
  70. A. Perrier, S. Aloise, M. Olivucci and D. Jacquemin, J. Phys. Chem. Lett., 2013, 4, 2190–2196 CrossRef CAS.
  71. D. P. Hoffman, S. R. Ellis and R. A. Mathies, J. Phys. Chem. A, 2014, 118, 4955–4965 CrossRef CAS.
  72. W. Liu, Y. Wang, L. Tang, B. G. Oscar, L. Zhu and C. Fang, Chem. Sci., 2016, 7, 5484–5494 RSC.
  73. C. Fang, R. R. Frontiera, R. Tran and R. A. Mathies, Nature, 2009, 462, 200–204 CrossRef CAS.

Footnote

Electronic supplementary information (ESI) available: The main structural parameters of the CHD molecule in ground and excited states; 2D wavelet maps from TRJ I, II and III. See DOI: 10.1039/c8cp04707c

This journal is © the Owner Societies 2019