Volume-conserving photoisomerization of a nonplanar GFP chromophore derivative: Nonadiabatic dynamics simulation

Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 214 (2019) 86–94

Contents lists available at ScienceDirect

Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy journal homepage: www.elsevier.com/locate/saa

Volume-conserving photoisomerization of a nonplanar GFP chromophore derivative: Nonadiabatic dynamics simulation Aihua Gao ⁎, Meishan Wang School of Physics and Optoelectronic Engineering, Ludong University, Yantai 264025, China

a r t i c l e

i n f o

Article history: Received 4 December 2018 Received in revised form 17 January 2019 Accepted 2 February 2019 Available online 5 February 2019 Keywords: GFP chromophore derivative Photoisomerization Surface hopping dynamics Reaction pathways

a b s t r a c t Nonadiabatic dynamics of a nonplanar green fluorescent protein (GFP) chromophore derivative was examined theoretically by trajectory surface hopping approach at the CASSCF level based on Zhu-Nakamura theory. The geometry optimizations show that there are four ground-state minima and four conical intersections between the ground (S0) and first excited (S1) states. Four S1-state minima were found at a perpendicularly twisted conformation around the imidazolinone-bridged bond. Upon excitation to S1 state, the main decay pathways of four isomers involve different S0/S1 potential energy surface conical intersections. The dominant excited-state relaxation mechanism of this GFP chromophore derivative is the twists of two bridging bonds and the methyl group in the bridge combined with the pyramidalization character of the central carbon atom. Further twists of two bridging bonds and the methyl group occur sequentially in the S0 state. It is worth to mention that the special volume-conserving motion of this molecule is attributed to twists of two bridging bonds in the same direction during the whole photoisomerization processes. The theoretical investigation presented herein provides important insights into the volume-conserving photoisomerization mechanisms of a nonplanar GFP chromophore derivative. We believe that the present work can benefit the design of the photochromic molecule devices in confined media. © 2019 Published by Elsevier B.V.

1. Introduction Since the green fluorescent protein (GFP) was first discovered in jellyfish Aquorea Victoria, it has been the subject of considerable studies and used as essential tools in molecular biology, biophysics, and biochemistry [1–11]. p-Hydroxybenzilideneimidazolinone (HBDI, Scheme 1) was widely used as the model GFP chromophore. In GFPs, the fluorescent on-state and nonfluorescent off-state was found to be anionic cis and neutral trans HBDI chromophore, respectively. Many studies have been carried out on the photophysical and photochemical properties of the simplest chromophore (HBDI) and its numerous derivatives [12–24]. It was reported that the radiationless decay from the first singlet excited state (S1) to ground state (S0) in HBDI occurs by twisting of two bridging bonds, which are connected to the imidazolinone and phenol rings (τI and τP twists, Scheme 1). A simultaneous rotation of the two bridging bonds (known as hula-twist) was considered as the most preferred pathway for the anionic cis form [14–17]. The intramolecular motions promoting radiationless decay in vacuo are very similar to those observed in solution [13]. Recently, Park et al. proposed two

⁎ Corresponding author. E-mail addresses: [email protected] (A. Gao), [email protected] (M. Wang).

https://doi.org/10.1016/j.saa.2019.02.001 1386-1425/© 2019 Published by Elsevier B.V.

main reaction pathways to explain the decay in cis HBDI [20]. The imidazolinone-bridged bond channel (τI) is suggested to be the major nonadiabatic decay pathway and the phenol-bridged bond channel (τP) is observed in about 30% of the trajectories in the on-the-fly XMS-CASPT2 (extended multistate multireference second-order perturbation theory) surface hopping simulations. For the neutral trans chromophore, the τI twist plays an important part in the internal conversion, while the τP twist is restricted in the ultrafast excited-state pathway [25]. The effort has been largely devoted to unravelling the photochemistry of GFP chromophore and its variants, in which the structure of the bridge, imidazolinone and phenol rings is usually planar. However, very little attention has been payed to the photoswitching mechanisms of the nonplanar chromophore in fluorescent proteins. Recently, Meech and Blancafort et al. have experimentally and theoretically investigated the photochemistry and photophysics of a novel sterically crowded nonplanar derivative of the GFP chromophore (brMe-HBDI, Scheme 1) [26]. The decay processes of brMe-HBDI are not dependent of solvent polarity and viscosity. This methylated derivative undergoes an excited-state photoisomerization with an ultrafast decay. After excited to S1 state, the system follows a barrierless, volume-conserving coordinate composed of torsion of two rings in the same direction and simultaneous pyramidalization of the methyne bridge carbon. In addition, a volume-conserving motion (or a pedalo-type motion) was found in an

A. Gao, M. Wang / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 214 (2019) 86–94

87

azodicarboxamide-based molecular switch in time-resolved infrared spectroscopy experiment and quantum-mechanical calculation [27]. It is noteworthy that the systems, such as brMe-HBDI or azodicarboxamide-based molecular switch, are able to influence the local environment with volume-conserving motions without the viscosity effects of media. Thus, these volume-conserving properties during the trans-cis or cis-trans photoisomerization facilitate the design of functional systems in confined media. The purpose of this article is to provide direct evidence of the volume-conserving motion in the photoisomerization of GFP chromophore derivative brMe-HBDI. From a theoretical point of view, it is necessary to perform full-dimensional nonadiabatic dynamics simulations to elucidate the photoinduced mechanisms of radiationless relaxation. In this work, complete active space self-consistent field (CASSCF) [28,29] and multiconfiguration second-order perturbation theory (CASPT2) [30] was used to obtain the features of S0 and S1 potential energy surfaces, which was followed by a surface hopping dynamics simulation. As a result, a clear picture of the photoisomerization dynamics of brMe-HBDI was obtained. The photoisomerization of brMe-HBDI follows a volume-conserving motion due to the torsion of two bridging bonds in the same direction.

Nonadiabatic dynamics of cis and trans brMe-HBDI photoisomerization was interpreted with the aid of trajectory surface hopping simulation based on the Zhu-Nakamura formula [32–35]. Initial conditions of the nuclear coordinates and velocities were given by a 200-ps molecular dynamics simulation of brMe-HBDI in the ground state at a temperature of 298 K. Classical trajectories showed nuclear motions and were propagated via a numerical integration of Newton's equations. The Newton's equations were integrated with the velocityVerlet method [36]. Time step was chosen to be 0.5 fs for adiabatic trajectory propagation and was set to 0.1 fs near conical intersection regions. All relevant electronic structures at each step were calculated on-the-fly using an external quantum chemistry program (Molpro 2010.1). The required energies and gradients were evaluated by two-state averaged CASSCF(6,5)/6-31G calculation, which was reported to be suitable for stable dynamics simulations for GFP chromophore [17,25]. The active space for CASSCF(6,5) was shown in Fig. S2. The simulation time for calculated trajectories was 1.5 ps and all trajectories dynamics calculations started from the S1 state. The nonadiabatic dynamics simulations in this work were carried out with NAIMD-DICP package, which has been successfully applied to understand the nonadiabatic decay processes in several systems [37–42].

2. Computational methods

3. Results and discussion

In this work, CASSCF calculations were carried out for brMe-HBDI to optimize the stationary points on S0 and S1 potential energy surfaces and conical intersections between two states. Ten electrons in eight orbitals were included in the active space (CASSCF(10,8), see Fig. S1) with the two electronic states at the same weights. The 6-31G(d) basis set was used. In addition, transition state structures obtained in CASSCF (10,8) calculation were characterized by analytical frequency calculations. In order to obtain more reliable energies, the single point energy correction calculations were performed using CASPT2 based on the CASSCF(10,8) optimized geometries. Relative and vertical energies of the CASSCF(10,8) optimized geometries were calculated at the CASPT2(10,8)/6-31G(d) level (CASPT2//CASSCF(10,8)/6-31G(d)). Additionally, an increased active space (14 electrons in 12 orbitals) was employed in the CASPT2 calculation (see Fig. S1). The excited-state potential energy curve calculations were performed at the CASPT2(14,12)/ 6-31G(d) level to investigate the excited-state pathway of cis and trans brMe-HBDI. All calculations were performed using the Molpro 2010.1 program [31].

3.1. Optimized geometric structures We performed an unconstrained optimization in the S0 and S1 states at CASSCF(6,5)/6-31G and CASSCF(10,8)/6-31G(d) levels. The dihedral angles N1-C2-C3-C4 (τI) and C2-C3-C4-C5 (τP) (see Scheme 1) represent the twists around imidazolinone and phenol bonds, respectively. The ground- and excited-state stable geometries of brMe-HBDI are shown in Fig. 1 and key geometrical parameters are presented in Table 1. The results of CASSCF(6,5)/6-31G level, which was used in the dynamics calculations, were compared with the results obtained in CASSCF(10,8)/6-31G(d). Here, the difference in the τP angle is around 20°–30° in the CASSCF calculations, while both levels of other geometrical parameter results agree well. In this work, we focused on the torsion motions of bridging bonds (τI and τP dihedral angles) and thus the difference in the τP angle is irrelevant for the subsequent dynamics simulations. There are two nonplanar stable structures for both cis and trans configurations in the S0 state. Compared with the planar form in HBDI, the interaction between the two rings and the methyl substituent in

Scheme 1. The cis-trans and trans-cis isomerization of GFP chromophore HBDI (upper) and its derivative brMe-HBDI (lower).

88

A. Gao, M. Wang / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 214 (2019) 86–94

Fig. 1. The ground- and excited-state minima, and conical intersections optimized at the CASSCF(10,8)/6-31G(d) level.

the bridge leads to the twists of imidazolinone and phenol bonds to form a stable minimum in brMe-HBDI. As a result, the chemical modification of the methyl group in the bridge removes the planar configuration that originally exists in HBDI. In particular, C3\\C4 has a more significant single bond character than C2\\C3 so that the phenolbridged bond twists more remarkably to drive the molecule toward a ground-state minimum (Table 1 and Fig. 1). Additionally, ground-state transition states were optimized to estimate the energy barriers

between two cis or trans isomers (Fig. S3 and Table 1). The barrier between two cis isomers is 0.002 eV and the barrier for two trans isomers is calculated to be 0.6 eV at the CASPT2 level. Thus, there is reversible thermal isomerization between two cis or trans isomers in the S0 state. For the four optimized S1-state minima, C2\\C3 bond elongates from 1.33–1.37 Å to 1.45 Å to make τI rotation more facile and the imidazolinone bond twists about 90°, giving rise to nearly a perpendicular conformation in the S1 state.

Table 1 Important geometrical parameters for optimized structures in CASSCF(6,5)/6-31G and CASSCF(10,8)/6-31G(d) calculation (bond lengths in angstroms and angles in degrees). Structure

Method

S0-cis1

CAS(6,5) CAS(10,8) CAS(6,5) CAS(10,8) CAS(6,5) CAS(10,8) CAS(6,5) CAS(10,8) CAS(6,5) CAS(10,8) CAS(6,5) CAS(10,8) CAS(6,5) CAS(10,8) CAS(6,5) CAS(10,8) CAS(10,8) CAS(10,8) CAS(10,8) CAS(10,8) CAS(10,8) CAS(10,8)

S0-cis2 S0-trans1 S0-trans2 S1-min1 S1-min2 S1-min3 S1-min4 TSab TScd CI1 CI2 CI3 CI4

C2-C3

C3-C4

N1-C2-C3

C2-C3-C4

C3-C4-C5

1.37 1.35 1.37 1.35 1.37 1.33 1.37 1.35 1.45 1.45 1.45 1.45 1.45 1.45 1.45 1.45 1.34 1.32 1.46 1.46 1.46 1.46

1.48 1.49 1.48 1.49 1.47 1.49 1.48 1.50 1.40 1.41 1.40 1.41 1.40 1.41 1.40 1.41 1.50 1.50 1.40 1.40 1.40 1.40

125.3 124.7 125.3 124.8 121.1 124.4 121.6 124.3 122.3 122.6 122.4 122.7 122.3 122.6 122.5 122.7 125.5 124.9 115.4 115.1 115.4 115.1

123.8 121.4 123.6 121.5 126.0 123.3 125.4 123.1 121.5 121.7 121.4 121.7 121.5 121.7 121.4 121.7 122.5 123.4 121.9 121.9 121.9 121.9

119.4 120.5 122.9 121.9 119.1 120.2 122.4 121.4 124.0 124.2 119.2 118.8 124.0 124.2 119.2 118.8 122.4 122.9 124.9 117.1 124.9 117.1

τI

τP

τpyr

2.3 1.9 2.6 1.9 −178.4 −178.9 −176.9 −178.7 93.3 95.1 93.1 95.0 −93.2 −95.0 −93.0 −95.0 −0.9 −179.2 108.7 109.1 −108.8 −109.1

−156.2 −131.7 26.4 48.5 −146.1 −114.5 39.4 71.1 179.4 178.8 −1.3 −1.6 −179.4 −178.8 1.3 1.6 −36.5 −22.1 −179.4 1.2 179.4 −1.3

0.03 0.04 0.08 0.00 0.07 0.06 0.04 0.09 0.12 0.24 0.12 0.25 0.12 0.24 0.13 0.25 0.08 0.03 0.37 0.38 0.37 0.38

A. Gao, M. Wang / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 214 (2019) 86–94

89

Table 2 Computed relative and vertical excitation energies (in eV) of ground- and excited-state minima. Energy

Method

Erea

CAS(6,5)b CAS(10,8)c CASPT2d CASPT2e CAS(6,5)b CAS(10,8)c CASPT2d CASPT2e

Everf

a b c d e f

S0-cis1

S0-cis2

S0-trans1

S0-trans2

S1-min1

S1-min2

S1-min3

S1-min4

0 0 0 0 5.14 5.65 4.38 4.57

0 0 0 0 5.23 5.64 4.53 4.56

0.10 0.0015 0.088 0.059 5.01 4.88 4.49 4.66

0.10 0.0014 0.080 0.090 5.13 5.71 4.52 4.91

2.93 3.57 2.69 2.77 0.26 0.75 0.87 1.04

3.00 3.61 2.71 2.77 0.33 0.78 0.90 1.00

2.93 3.57 2.69 2.77 0.26 0.75 0.87 1.04

3.00 3.61 2.71 2.77 0.34 0.78 0.90 1.00

Energies relative to S0-cis1. 6-31G. 6-31G(d). CASPT2(10,8)/6-31G(d). CASPT2(14,12)/6-31G(d). Vertical energies.

Minimum energy conical intersections (CIs) between potential energy surfaces play a crucial role in the internal conversion processes [43,44]. Four S0/S1 CIs for brMe-HBDI were located at the CASSCF (10,8)/6-31G(d) level (Fig. 1). The key geometrical parameters are collected in Table 1. Structurally, the optimized CI1 (CI2) is a chiral isomer of CI3 (CI4). The dihedral angles τI are around ±109°, indicating that the molecule needs a further imidazolinone bond twist (about 15°) from a S1-state minimum to a CI in the S1 state for two cis isomers. The CI structures, CI1 and CI2, are consistent with previous studies for brMe-HBDI [26]. For the CIs of original neutral HBDI chromophore, the dihedral angles, N1-C2-C3-C4 and C2-C3-C4-C5, are around ±85° and ±18°, respectively [25]. Based on the optimized geometries, both the S1minima and CIs are characterized by the τI and τP rotation in the excited state. In our work, pyramidalization at the central carbon atom (C3) is described by τpyr (τpyr = (360° − S)1/2, S is the sum of C2\\C3\\C6, C2\\C3\\C4, and C4\\C3\\C6 bond angle values) [45]. According to the values of τpyr for the optimized geometries in Table 1, the atoms, C2, C3, C4, and C6, are nearly in a planar structure for S0-minima. The pyramidalization characteristic of C3 is observed for S1-minima and CIs. Table 1 shows that the geometries of the optimized S0/S1 CIs are characterized not only by the τI and τP twisting but also by the pyramidalization of the C3 atom. Comparing the results of theoretical work on CIs in HBDI chromophore anion (3.0–3.8 degree1/2) [17], the pyramidalization characteristic is weak for the CIs of brMe-HBDI (0.37–0.38 degree1/2). Tables 2 and 3 list the relative (referred to S0-cis1 and S0-cis2) and vertical excitation energies at CASSCF(6,5)/6-31G, CASSCF(10,8)/631G(d), CASPT2(10,8)/6-31G(d), and CASPT2(14,12)/6-31G(d) levels. The vertical excitation energies for two cis isomers are evaluated to be 4.38–4.57 eV at the CASPT2 level. Experimental measurements show that irradiation into S1 with 365 ± 20 nm light converts cis to trans

Table 3 Computed relative and vertical excitation energies (in eV) of conical intersections. Energy

Method

CI1

CI2

CI3

CI4

Erea

CAS(10,8)b CASPT2c CASPT2d CAS(10,8)b CASPT2c CASPT2d

3.84 2.81 2.96 0.0003 0.27 0.34

3.90 2.91 2.99 0.002 0.23 0.38

3.84 2.81 2.94 0.002 0.27 0.35

3.90 2.91 2.89 0.0004 0.23 0.24

Evere

a b c d e

Energies relative to S0-cis1. 6-31G(d). CASPT2(10,8)/6-31G(d). CASPT2(14,12)/6-31G(d). Vertical energies.

isomers, and cis and trans isomers have similar spectra [26]. It was reported that the results of CASPT2 levels are accurate within 0.1–0.2 eV [46]. In Table 2, both cis and trans isomers are very close energetically at CASPT2 levels, consistent with the results of MP2 calculation [26]. The CASPT2 calculations indicate that the energetic level of four S1state minima is almost equal (2.77 eV at CASPT2(14,12)/6-31G (d) level), in agreement with the experimentally measured emission maximum wavelength around 427 nm (2.90 eV) [26]. The fluorescence emission was reported to be obtained by inversion of the central bonds [26]. However, based on the perpendicular conformation of S1-state minima (Table 1), the emission in this molecule is related to relaxation processes involving not only the stretching but also the torsion of two central bonds. The CIs lie 0.10–0.22 eV above the S1-state minima. The energy gaps of CIs were calculated b0.38 eV at CASPT2 level to confirm our location of S0/S1 conical intersections. Moreover, we calculated the excited-state energy profiles between the Franck-Condon region of S0-cis1 and S1-min1 at CASSCF(6,5)/6-31G, CASSCF(10,8)/6-31G(d), and CASPT2(10,8)/6-31G(d) levels of theory (see Fig. S4). The results provide the available adequacy of CASSCF calculations. 3.2. Nonadiabatic dynamics of cis and trans brMe-HBDI All cis and trans isomers were considered in the nonadiabatic dynamics simulations and a total of 40 trajectories were run for the photoisomerization reactions of brMe-HBDI. The trajectories are initiated at geometries in the vicinity of the respective Franck-Condon region of S0-cis1, S0-cis2, S0-trans1, and S0-trans2 in the S1 state. Among the 20 trajectories for cis isomers, isomerization to trans isomer occurs in 10 trajectories. In 20 trajectories simulated for trans isomers, isomerization to cis form occurs in 14 trajectories. Therefore, the isomerization quantum yields of 0.50 and 0.70 for cis and trans isomers are obtained, respectively. These results are not consistent with the experimental values (4 ± 2% and 25 ± 10% for cis and trans isomers, respectively) [26]. The theoretical larger quantum yields may be caused by the small number of trajectories or neglection of solvent effects in the dynamics simulations. Further experimental or theoretical studies are necessary for photoconversion quantum yields of this system. In current work, we just focus on the photoinduced isomerization mechanisms of brMe-HBDI. According to our optimization, the τI twist around imidazolinone bond and τP around phenol bond may be involved in the radiationless relaxation. Thus, the dihedral angles N1-C2-C3-C4 (τI) and C2-C3-C4C5 (τP) of the hopping geometries at which surface hops occur are summarized for all trajectories in Fig. 2. Obviously, the hopping points of all trajectories are around the optimized CI regions, showing the importance of the CIs in the decay process for this molecule. The surface hops in the dynamics of two cis isomers take place mainly around CI1

90

A. Gao, M. Wang / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 214 (2019) 86–94

Fig. 2. Geometrical distributions of the dihedral angles C2-C3-C4-C5 (τP) and N1-C2-C3-C4 (τI) at the S1/S0 hopping points. Black circles and dots denote hopping points for S0-cis1 and S0-cis2, respectively. Green circles and dots show hopping points for S0-trans1 and S0-trans2, respectively. Red triangles represent the four optimized CIs.

and CI2, while the hops occur predominantly around CI3 and CI4 for two trans isomers (see Table 1). To gain more detailed structural insight into the photoisomerization reactions, we analyze the variations of key geometric parameters during the whole dynamics simulation. Typical trajectories for isomerization in S0-cis1 and S0-cis2 are shown in Fig. 3(a,b) and (c,d), respectively. In the trajectory for S0-cis1, the system evolves on the S1 potential energy surface and the decay to S0 occurs at 250 fs. During the initial stage (b15 fs), the C2\\C3 bond shows an obvious elongation from 1.35 Å to 1.55 Å to increase single bond character, enhancing the τI twist in the

S1 state. Meanwhile, the C3\\C4 bridging bond becomes shorter but still has a single bond character. In Fig. 3(b), it can be concluded that both the τI and τP twists are involved in the excited-state pathway. At the hopping point, τI and τP are nearly 85° and 170°, respectively. Here, a twisted conformation is generated around the C2-C3 bond corresponding to a conical intersection (CI1, see Table 1). Then, the molecule decays to the ground state. In the S0 state, further τI and τP twists around two bridging bonds make the system to achieve a trans isomer region (S0-trans2) at about 400 fs. Subsequently, further τP twist leads the system to cross a small potential energy barrier between two trans isomers, approaching to the S0-trans1 region at about 1100 fs. Additionally, the C2-C3-C6-H7 dihedral angle (τM) was used to describe the behavior of the methyl group linked to the central carbon C3. In Fig. 3(b), a swinging motion of the methyl group is observed during the whole simulation, which is not too surprising since the molecule has a sterically crowded structure and the twist of the methyl group is necessary for the reorientation of the imidazolinone and phenol rings. This phenomenon was predicted by calculations combining TD-CAMB3LYP and CASPT2 method [26]. In the photoisomerization process of S0-cis2 (see Fig. 3(c,d)), two bridging bonds exhibit the same behaviors as those in Fig. 3(a,b). The hopping event occurs in the vicinity of CI2. The system reaches a trans isomer region (S0-trans1) at about 500 fs. The τM twist is also found in this trajectory. In a word, the photoisomerization reactions for two cis isomers are characterized not only by the twisting of two bridging bonds (τI and τP) but also by the twisting of the methyl group (τM). According to the optimized geometries shown in Table 1, upon excitation to the S1 state, the system could relax from the Franck-Condon region of S0-cis1 or S0-cis2 to a S1-state minimum, S1-min1 or S1-min2, and then decay to the ground state via a CI. Fig. 4 shows typical trajectories for photoisomerization processes of S0-trans1 (top two panels) and S0-trans2 (bottom two panels). After excitation to the S1 state, the behaviors of two bridging bonds are found to be similar with those observed in the photoisomerization reactions of two cis isomers. In Fig. 4(b), the molecule reaches a conical intersection region (CI3, see Table 1) by the twists of two bridging

Fig. 3. Variations in the bond lengths and dihedral angles in the typical trajectories (a,b) for S0-cis1 and (c,d) for S0-cis2. The vertical red line marks the S1/S0 hopping event.

A. Gao, M. Wang / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 214 (2019) 86–94

91

Fig. 4. Variations in the bond lengths and dihedral angles in the typical trajectories (a,b) for S0-trans1 and (c,d) for S0-trans2. The vertical red line marks the S1/S0 hopping event.

bonds in the S1 state. Near the CI3 region at 540 fs, the τI and τP dihedral angles are nearly −90° and 145°, respectively, and the nonadiabatic decay to the S0 state takes place. In the S0 state, the system reaches a cis isomer region (S0-cis2). Significantly, the twists of two bridging bonds are also observed in the whole trajectory in Fig. 4(d). The system decays to the S0 state through CI4 and achieves the S0-cis1 region at around 600 fs. It is worth to mention that the τM twist is also necessary in the isomerization processes of two trans isomers. In all simulated trajectories, the photoisomerization of trans isomers excited to S1 needs longer time scale than that of cis isomers. Specifically, the photoisomerization starting from a trans conformer does not involve S1-state minima in the excited-state pathway. Fig. 5 shows the pyramidalization feature of the C3 atom at hopping events for all trajectories. The time evolution of pyramidalization angles for trajectories in Figs. 3 and 4 is presented in Fig. S5. It can be

demonstrated that the hopping events take place predominantly with the pyramidalization angles ranging form 0.5 to 1.7 degree1/2. Thus, the twists of two bridging bonds and the methyl group are coupled by the pyramidalization of the C3 atom. Han et al. reported a larger pyramidalization angle at the hopping events in the nonadiabatic dynamics of HBDI anion [17]. Our conclusion is that, to a certain extent, the interaction between the methyl group and two rings could hinder the pyramidalization of central carbon atom in brMe-HBDI. The large scale bond torsion during the isomerization may be opposed by solvent environments. Nevertheless, Meech and Blancafort et al. found that this methylated derivative, brMe-HBDI, shows an exceptionally fast excited state decay which is essentially independent of solvent viscosity [26]. In this work, we evaluated the structural changes by calculating the molecular length in the photoisomerization reactions of this system. Fig. 6 shows the time evolution of molecular length in the isomerization processes in four typical trajectories. The cis, trans isomers, and CIs are predicted to have similar length (10.6–11.0 Å). As shown in four panels, it is interesting to note that the isomerization occurs with small changes in the length (b0.9 Å), indicating that the isomerization reaction of the target molecule turns out to follow a special volume-conserving motion. To explore this volume-conserving property, we analyze the directions of two rings rotations. The initial structures of typical trajectories and the respective structures at the hopping events are presented in Fig. 7. As shown, the two bridging bonds are rotating in the same direction upon the chromophore excited to the S1 state. These results and the time evolution of the dihedral angles (see Figs. 3 and 4) reveal that two bridging bonds rotate in the same direction during the whole photoisomerization processes of brMeHBDI. Hence this twist motion leads to the volume-conserving coordinate and the lack of a solvent viscosity effect for the internal conversion. It is in consistent with the static quantum calculation in Ref. [26]. 3.3. Excited-state pathways of cis and trans brMe-HBDI

Fig. 5. Geometrical distribution of the pyramidalization angles at the S1/S0 hopping points.

The decay as a function of time is shown in Fig. 8. In the end, all trajectories decay to the S0 state. In the decay processes for cis isomers,

92

A. Gao, M. Wang / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 214 (2019) 86–94

Fig. 6. Time evolution of the length of the molecule during the photoisomerization for (a) S0-cis1, (b) S0-cis2, (c) S0-trans1, (d) S0-trans2.

most of the S1/S0 hopping events take place between 170 and 400 fs, which is in good agreement with the experimental results of exponentially ultrafast decay (faster than 400 fs) [26]. Moreover, it can be seen that the S1-state lifetime for the photoinduced isomerization reaction of trans isomers is a little longer than that of the cis isomers. We calculated the excited-state energy profiles along the constructed linearly interpolated internal coordinate pathway at the CASPT(14,12)/6-31G (d) level (see Fig. S6). For S0-cis1, the pathway is barrierless toward a S1-state minimum point. After photoexcitation, S0-trans1 could relax

from the Franck-Condon point to the CI3 region on the potential energy surface in a barrierless manner. Specifically, the longer decay time for trans isomers can be attributed to a flat region in the initial relaxation pathway. 4. Conclusion By using quantum calculations and surface hopping molecular dynamics simulations, we probe the photoisomerization mechanisms of

Fig. 7. Initial structures of typical trajectories (in white) and the respective structures at the hopping events in the nonadiabatic dynamics simulations of two cis (upper) and trans isomers (lower).

A. Gao, M. Wang / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 214 (2019) 86–94

93

Acknowledgments This work was supported by National Natural Science Foundation of China (Grants 11704171) and the Open Fund of the State Key Laboratory of Molecular Reaction Dynamics in DICP, CAS. Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi. org/10.1016/j.saa.2019.02.001.

References

Fig. 8. The decay of the S1 population as a function of time derived from nonadiabatic dynamics simulations.

brMe-HBDI, a nonplanar GFP chromophore derivative. In the optimization, two cis and two trans isomers are found in the S0 state and there are four S0/S1 CIs between two involved states. Moreover, four optimized S1-state minima are predicted to have closely perpendicular conformations. In terms of quantum calculations and dynamics simulations, we summarized the photoisomerization mechanisms of four groundstate minima in brMe-HBDI in Fig. 9. Upon photoexcitation to the S1 state, the radiationless relaxation of four isomers mainly proceeds via different conical intersections. The twists of two bridging bonds (τI and τP) and the methyl group (τM) combined with the pyramidalization of C3 atom (τpyr) are suggested to be the dominant mechanism in the photoisomerization processes. The target molecule needs further rotation of τI, τP, τM to reach the product region in the S0 state. The twists of two bridging bonds are in the same direction, leading to a special volume-conserving motion during the whole photoisomerization processes. Additionally, the computed decay lifetime of cis isomers is estimated about 400 fs, which is in good agreement with the experimentally measured value. However, a longer lifetime is observed for trans isomers due to a flat region on S1-state potential energy surface. Our work provides the detailed mechanisms of a special volumeconserving motion in a nonplanar GFP chromophore derivative and we hope that this work could be useful for the design of photochromic molecule devices in confined environments.

[1] F.H. Johnson, O. Shimomura, Y. Saiga, L.C. Gershman, G.T. Reynolds, J.R. Waters, Quantum efficiency of cypridina luminescence, with a note on that of Aequorea, J. Cell. Comp. Physiol. 60 (1962) 85. [2] R.Y. Tsien, The green fluorescent protein, Annu. Rev. Biochem. 67 (1998) 509. [3] A. Miyawaki, Innovations in the imaging of brain functions using fluorescent proteins, Neuron 48 (2005) 189. [4] S. Rafiq, B.K. Rajbongshi, N.N. Nair, P. Sen, G. Ramanathan, Excited state relaxation dynamics of model green fluorescent protein chromophore analogs: evidence for cis-trans isomerism, J. Phys. Chem. A 115 (2011) 13733. [5] J.S. Paige, K.Y. Wu, S.R. Jaffrey, RNA mimics of green fluorescent protein, Science 333 (2011) 642. [6] J.S. Donner, S.A. Thompson, M.P. Kreuzer, G. Baffou, R. Quidant, Mapping intracellular temperature using green fluorescent protein, Nano Lett. 12 (2012) 2107. [7] K. Addison, I.A. Heisler, J. Conyard, T. Dixon, P.C.B. Page, S.R. Meech, Ultrafast excited state dynamics of the green fluorescent protein chromophore and its kindling fluorescent protein analogue, Faraday Discuss. 163 (2013) 277. [8] Z.G. Chiragwandi, K. Gillespie, Q.X. Zhao, M. Willander, Ultraviolet driven negative current and rectifier effect in self-assembled green fluorescent protein device, Appl. Phys. Lett. 89 (2006), 162909. . [9] V. Adam, H. Mizuno, A. Grichine, J.I. Hotta, Y. Yamagata, B. Moeyaert, G.U. Nienhaus, A. Miyawaki, D. Bourgeois, J. Hofkens, Data storage based on photochromic and photoconvertible fluorescent proteins, J. Biotechnol. 149 (2010) 289. [10] D.M. Shcherbakova, P. Sengupta, J. Lippincott-Schwartz, V.V. Verkhusha, Photocontrollable fluorescent proteins for superresolution imaging, Annu. Rev. Biophys. 43 (2014) 303. [11] D.M. Chudakov, M.V. Matz, S. Lukyanov, K.A. Lukyanov, Fluorescent proteins and their applications in imaging living cells and tissues, Physiol. Rev. 90 (2010) 1103. [12] R.S.H. Liu, G.S. Hammond, The case of medium-dependent dual mechanisms for photoisomerization: one-bond-flip and hula-twist, Proc. Natl. Acad. Sci. U. S. A. 97 (2000) 11153. [13] C.R.S. Mooney, D.A. Horke, A.S. Chatterley, A. Simperler, H.H. Fielding, J.R.R. Verlet, Taking the green fluorescence out of the protein: dynamics of the isolated GFP chromophore anion, Chem. Sci. 4 (2013) 921. [14] K.L. Litvinenko, N.M. Webber, S.R. Meech, Internal conversion in the chromophore of the green fluorescent protein: temperature dependence and isoviscosity analysis, J. Phys. Chem. A 107 (2003) 2616. [15] D. Mandal, T. Tahara, S.R. Meech, Excited-state dynamics in the green fluorescent protein chromophore, J. Phys. Chem. B 108 (2004) 1102.

Fig. 9. Proposed photoisomerization mechanisms of brMe-HBDI. The dashed and solid arrows denote the main excited- and ground-state reaction channels, respectively. Black and red arrows show the channels for two cis isomers. The pathways of two trans isomers are marked by blue and green arrows.

94

A. Gao, M. Wang / Spectrochimica Acta Part A: Molecular and Biomolecular Spectroscopy 214 (2019) 86–94

[16] M.E. Martin, F. Negri, M. Olivucci, Origin, nature, and fate of the fluorescent state of the green fluorescent protein chromophore at the CASPT2//CASSCF resolution, J. Am. Chem. Soc. 126 (2004) 5452. [17] L. Zhao, P. Zhou, B. Li, A. Gao, K. Han, Non-adiabatic dynamics of isolated green fluorescent protein chromophore anion, J. Chem. Phys. 141 (2014), 235101. . [18] S. Olsen, K. Lamothe, T.J. Martínez, Protonic gating of excited-state twisting and charge localization in GFP chromophores: a mechanistic hypothesis for reversible photoswitching, J. Am. Chem. Soc. 132 (2010) 1192. [19] A. Toniolo, S. Olsen, L. Manohar, T.J. Martínez, Conical intersection dynamics in solution: the chromophore of green fluorescent protein, Faraday Discuss. 127 (2004) 149. [20] J.W. Park, T. Shiozaki, On-the-fly CASPT2 surface-hopping dynamics, J. Chem. Theory Comput. 13 (2017) 3676. [21] L.M. Tolbert, A. Baldridge, J. Kowalik, K.M. Solntsev, Collapse and recovery of green fluorescent protein chromophore emission through topological effects, Acc. Chem. Res. 45 (2012) 171. [22] K.B. Bravaya, B.L. Grigorenko, A.V. Nemukhin, A.I. Krylov, Quantum chemistry behind bioimaging: insights from ab initio studies of fluorescent proteins and their chromophores, Acc. Chem. Res. 45 (2012) 265. [23] P.J. Tonge, S.R. Meech, Excited state dynamics in the green fluorescent protein, J. Photochem. Photobiol., A 205 (1) (2009). [24] Y.H. Hsu, Y.A. Chen, H.W. Tseng, Z. Zhang, J.Y. Shen, W.T. Chuang, T.C. Lin, C.S. Lee, W.Y. Hung, B.C. Hong, S.H. Liu, P.T. Chou, Locked ortho- and para-core chromophores of green fluorescent protein; dramatic emission enhancement via structural constraint, J. Am. Chem. Soc. 136 (2014) 11805. [25] A. Gao, M. Wang, J. Ding, Ultrafast trans-cis photoisomerization of the neutral chromophore in green fluorescent proteins: surface-hopping dynamics simulation, J. Chem. Phys. 149 (2018), 074304. . [26] J. Conyard, I.A. Heisler, Y. Chan, P.C.B. Page, S.R. Meech, L. Blancafort, A new twist in the photophysics of the GFP chromophore: a volume-conserving molecular torsion couple, Chem. Sci. 9 (2018) 1803. [27] S. Amirjalayer, A. Martinez-Cuezva, J. Berna, S. Woutersen, W.J. Buma, Photoinduced pedalo-type motion in an azodicarboxamide-based molecular switch, Angew. Chem. 130 (2018) 1810. [28] P. Siegbahn, A. Heiberg, B. Roos, B. Levy, A comparison of the super-CI and the Newton-Raphson scheme in the complete active space SCF method, Phys. Scr. 21 (1980) 323. [29] B.O. Roos, The complete active space SCF method in a fock-matrix-based super-CI formulation, Int. J. Quantum Chem. 18 (1980) 175. [30] K. Andersson, P. Malmqvist, B.O. Roos, Second-order perturbation theory with a complete active space self-consistent field reference function, J. Chem. Phys. 96 (1992) 1218.

[31] H.J. Werner, P.J. Knowles, G. Knizia, F.R. Manby, M. Schütz, et al., molpro, version 2010.1, a package of ab initio programs, (see) http://www.molpro.net 2010. [32] C. Zhu, H. Nakamura, The two-state linear curve crossing problems revisited. III. Analytical approximations for stokes constant and scattering matrix: nonadiabatic tunneling case, J. Chem. Phys. 98 (1993) 6208. [33] C. Zhu, H. Nakamura, Theory of nonadiabatic transition for general two-state curve crossing problems. I. Nonadiabatic tunneling case, J. Chem. Phys. 101 (1994) 10630. [34] C. Zhu, K. Nobusada, H. Nakamura, New implementation of the trajectory surface hopping method with use of the Zhu-Nakamura theory, J. Chem. Phys. 115 (2001) 3031. [35] P. Oloyede, G. Mil'nikov, H. Nakamura, Generalized trajectory surface hopping method based on the Zhu-Nakamura theory, J. Chem. Phys. 124 (2006), 144110. . [36] W.C. Swope, H.C. Andersen, P.H. Berens, K.R. Wilson, A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: application to small water clusters, J. Chem. Phys. 76 (1982) 637. [37] A.H. Gao, B. Li, P.Y. Zhang, K.L. Han, Nonadiabatic ab initio molecular dynamics of photoisomerization in bridged azobenzene, J. Chem. Phys. 137 (2012), 204305. . [38] A. Gao, M. Wang, Nonadiabatic ab initio molecular dynamics study of photoisomerization in N-salicilydenemethylfurylamine(SMFA), J. Chem. Phys. 146 (2017), 124312. . [39] A.H. Gao, B. Li, P.Y. Zhang, J.Y. Liu, Photochemical dynamics simulations for trans-cis photoisomerizations of azobenzene and bridged azobenzene, Comput. Theor. Chem. 1031 (2014) 13. [40] L. Zhao, J. Liu, P. Zhou, On the mechanism of non-radiative decay of blue fluorescent protein chromophore: new insight from the excited-state molecular dynamics simulations and potential energy calculations, Spectrochim. Acta A 186 (2017) 52. [41] L. Zhao, J. Liu, P. Zhou, New insight into the photoisomerization process of the salicylidene methylamine under vacuum, J. Phys. Chem. A 120 (2016) 7419. [42] L. Zhao, J. Liu, P. Zhou, Effect of methylation on the photodynamical behavior of arylazoimidazoles: new insight from theoretical ab initio potential energy calculations and molecular dynamics simulations, J. Phys. Chem. A 121 (2017) 141. [43] F. Bernardi, M. Olivucci, M.A. Robb, Potential energy surface crossings in organic photochemistry, Chem. Soc. Rev. 25 (1996) 321. [44] B.G. Levine, T.J. Martínez, Isomerization through conical intersections, Annu. Rev. Phys. Chem. 58 (2007) 613. [45] S.P. Gavrish, Approximate expressions and the relationship between pyramidalization (out-of-plane deformation) characteristics of trigonal centers, J. Comput. Chem. 33 (2012) 2173. [46] M. Merchán, L. Serrano-Andrés, Ab initio methods for excited states, in: M. Olivucci (Ed.), Comput. Photochem, Elsevier B.V., 2005.