Theoretical study of cyclohexadiene/hexatriene photochemical interconversion using spin-Flip time-Dependent density functional theory

Abstract

The photochemical interconversion between 1,3-cyclohexadiene (CHD) and all-cis-hexatriene (cZc-HT) is reinvestigated using spin-flip time-dependent density functional theory in combination with various hybrid functionals, BHHLYP functional showing the best performance. The critical geometries of the ground, S${}_0$, and the first two excited-state, S${}_1$ and S${}_2$, potential energy surfaces, such as, various minima, transition state, minimum-energy crossing points between S${}_2$/S${}_1$ and S${}_1$/S${}_0$ show an excellent agreement with those obtained by multireference wave function methods. Our results show how a low-cost method based on DFT can successfully describe and characterise the most important geometries on the potential energy surfaces along the ring-opening/closure reaction coordinate involved in the CHD to cZc-HT photoconversion.

GRAPHICAL ABSTRACT

Keywords:

• Photochemistry
• electrocyclic reactions
• spin-Flip
• conical intersection
• minimum-energy crossing point

1. Introduction

The photochemical interconversion between 1,3-cyclohexadiene (CHD) and all-cis-hexatriene (cZc-HT), that plays a crucial role in the photobiological synthesis of vitamin D${}_3$ [1,2], is one of the most widely studied photochemical electrocyclic reactions [3,4]. The nature of this photoconversion has been the subject of a number of both experimental and theoretical investigations [3–9], and it is often considered as simplified model for the investigation of the photophysical and photochemical properties of important macromolecules, that contain such (4n+2) photo electrocyclic reaction as their central building block. These photoactive molecules are relevant for modern technologies, for example, application of molecular photoswitches [10,11], photomolecular motors [11], nanomachines [12] in photomedicine [13] and molecular electronic devices [14,15].

The present status of experimental and theoretical studies on the photochemistry of CHD has been already reviewed by several authors [5,6,9,16–23] and many aspects of the CHD to cZc-HT photochemical interconversion have been clarified. Much of these works have focused on elucidating the ultrafast photoinduced dynamics of CHD through critical points on the potential energy surfaces, such as conical intersections (CIs). Even though it appears relatively simple, this is a very complex subject that led to numerous controversies in the field [9,19,22,24,25] and there exist contradictory experimental and theoretical findings [5,16,17,19,22]. However, it is now widely accepted that CHD to cZc-HT photochemical interconversion occurs via three consecutive steps [16,18,19,21,22] (see Figure 1); (i) Photoexcitation of the CHD from the ground state S${}_0$ (belongs to the irreducible representation 1${}^1$A under C${}_2$ symmetry) to the first excited state S${}_1$ (belongs to the irreducible representation 1${}^1$B under C${}_2$ symmetry) at the Franck-Condon (FC) region. (ii) Rapid decay from the excited state 1${}^1$B to the “dark” excited state 2${}^1$A. (iii) Decay from 2${}^1$A to either the ground state of the cZc-HT or CHD. To the best of our knowledge, the photoreversion of cZc-HT to CHD has not been observed, since cZc-HT structure does not remain in the thermal equilibrium, and subsequently derive into the predominant trans-HT structures via the thermal reactions [6].

Figure 1. Schematic representation of the CHD to cZc-HT photochemical interconversion. The dashed lines represent the crossing between excited states S${}_2$(2${}^1$A) and S${}_1$(1${}^1$B), and S${}_1$(2${}^1$A) and the ground state S${}_0$(1${}^1$A). FC refers to the Franck-Condon region.

The progression of the CHD electrocyclic ring opening described by this scheme is mediated mainly by two CIs, S${}_2$/S${}_1$ and S${}_1$/S${}_0$. Theoretical studies have been performed using multireference wave function methods [16,18,22], density functional theory (DFT) and time-dependent density functional theory (TDDFT) [19,21] in order to locate and characterise these important geometries. The complete active space self-consistent field method (CASSCF) is one of the multireference methods typically used to study the ring-opening/closure photoreaction between CHD and cZc-HT [3,7–9,26]. Usually, the active space employed with CASSCF consists of 6 electrons and 6 orbitals (4$\pi ,{\pi }^{\ast }$ and 2$\sigma ,{\sigma }^{\ast }$) for CHD and 6 electrons and 6 orbitals (6$\pi ,{\pi }^{\ast }$) for cZc-HT [9]. Although this method has satisfactorily described the reaction path between the states 2A and 1A [7,9,26], it fails to reproduce the correct order of the states 1B and 2A in the FC region [5,7–9]. This fact has limited the CASSCF method to describe this reaction process approximately after the crossing between S${}_2$ and S${}_1$ [8,9]. Santolini et al. [27] have obtained the correct order of the states in the FC region for the butadiene-like molecules by significantly increasing the size of the restricted active space (RAS). Although this demonstrates that the correct order of the states can be qualitatively achieved, the process of obtaining them for more complicated systems can be hard, and demands an increase in computational complexity [22,27]. The second-order perturbation theory for CASSCF method (CASPT2) is another widely used multireference method that has shown very good results [5,18,22]. Unlike CASSCF, this method computes the correct order of the states in the FC region. Mori and Kato [18] show an interesting comparison between extensions based on CASSCF and CASPT2. They concluded that the dynamic electron correlation is an important factor for the correct description of the energies and geometries [18]. Recently, Poliak et al. [22] used an extended multistate CASPT2, XMS-CASPT2, that offers a good balance of the dynamic electron correlation, for surface hopping simulations and provided a detailed description of the photodeactivation mechanism.

The DFT method, in contrast to multireference methods, has the ability to describe the dynamic electron correlation [28]. A previous work based on this method could capture the right order of states in the FC region [29], which were consistent with the results using CASPT2 [7]. However, DFT, and therefore TDDFT can not describe the nondynamic correlation [28]. Additionally, TDDFT shows serious problems to describe CIs, especially the CIs that involve the S${}_0$ [30,31]. Recently, Filatov et al. [21], used state-interaction, state-averaged spin-restricted ensemble-referenced Kohn-Sham (SSR) method, and were able to provide a correct description of the energies, geometries of the S${}_0$, and the CI${}_{S_1/S_0}$. However, the CI${}_{S_2/S_1}$ was not described since their studies were confined to the dynamic of CHD ring opening on the S${}_1$ state [21]. An extension to TDDFT proposed by Shao et al. [31], called spin-flip time-dependent density functional theory (SF-TDDFT), developed to describe diradicals with strong nondynamical correlation, has shown a great performance for describing CIs in ethylene [32], uracil [33] and cis-stilbene [34]. Within this approach, a high-spin triplet state is chosen as the initial reference state allowing that the S${}_0$ and the single excited states to be treated at the same footing. SF-TDDFT recovers both nondynamical and dynamical correlation from SF and DFT, respectively [31,35]. These features have made SF-TDDFT a great alternative to address photochemical processes. Here, we use SF-TDDFT to study and determine the important geometries at the potential energy profile of CHD to cZc-HT photochemical interconversion, with a special focus on locating and characterising the S${}_2$/S${}_1$ and S${}_1$/S${}_0$ CIs.

This paper is organised as follows: Section 2 briefly describes the underlying theoretical methods used in this work. Computational details are summerised in Section 3. Our results and discussions are presented in Section 4 and finally our concluding remarks are given in Section 5.

2. Theoretical methods

2.1. Orbital and state correlation diagrams

The CHD to cZc-HT photochemical interconversion can be depicted qualitatively using the conservation of orbital symmetry, i.e. following the Woodward-Hoffmann rules [36] for a pericyclic reaction. Similar to the description made for the ring closure of 1,3-butadiene in reference [37] using the Woodward-Hoffmann rules, we shall illustrate the photochemical interconversion between CHD and cZc-HT. Figure 2 describes an orbital correlation diagram when the absorption of a photon by the CHD leads to the excitation of an electron from the bonding orbital 2π to the antibonding orbital 3π. The conrotatory (4n+2 electron in π system) path of the excited CHD leads to the excited cZc-HT keeping the same orbital symmetry [37,38]. However, the suggestion of an intermediate excited state cZc-HT generates a problem since it is known that CHD, after photoexcitation, return to either its ground state or to the ground state of cZc-HT, i.e. without the formation of an intermediate excited state cZc-HT [4]. This discrepancy can be resolved if one considers the overall states of CHD and cZc-HT via state correlation diagram presented in Figure 3. For doing so, we used the direct product between the symmetry specie of the individual occupied orbitals satisfying the following rules: S × S = S; S × A = A and A × A = S, where S (with character 1) means symmetric with respect to a symmetry element and A (with character -1) means antisymmetric with respect to a symmetry element. All doubly occupied orbitals are always totally symmetric. Because the conrotatory path preserves the rotation around the C${}_2$ axis, we can classify the relevant states in terms of S and A with respect to this axis [37,38]. Since CHD (1σ,1π,2π) is correlated with cZc-HT (1π,2π,4π) (see Figure 2), it follows that the CHD (1${\sigma }^2$ 1${\pi }^2$ 2${\pi }^2$, ${}^1$S) is correlated with cZc-HT (1${\pi }^2$ 2${\pi }^2$ 4${\pi }^2$, ${}^1$S), where ${}^1$S denotes the result of the direct product between the symmetry S × S = S and the superscript 1 means that the electrons are in a singlet configuration (see Figure 3).

Figure 2. Orbital correlation diagram for the photochemical interconversion of the CHD to cZc-HT. The orbitals of CHD and cZc-HT are ordered in energy and labelled as symmetric (S) and antisymmetric (A) under C${}_2$ symmetry. Dashed lines relate empty orbitals with the same symmetry and solid lines relate occupied orbitals with the same symmetry.

Figure 3. Scheme of the state correlation diagram for the photocyclisation of CHD to cZc-HT. Dashed lines connect molecular orbitals with the same symmetry and same multiplicity, and solid blue lines connect molecular orbitals satisfying the non-crossing rule.

Following the same protocol in reference [37] for the ring closure of 1,3-butadiene example, the dashed lines in Figure 3 stand for the connections between the states with the same symmetry and multiplicity. Note that it is not allowed to have crossings between the energy levels with the same symmetry (non-crossing rule) [37,39,40], therefore the blue lines are used to remark what happens. However, the non-crossing rule is only valid for diatomic molecules. Thus, in polyatomic molecules, this rule loses its validity and two electronic states with the same symmetry and multiplicity may cross at a conical intersection [26,41,42]. Therefore, the dashed lines show the behaviour expected.

In the literature [3,17], the molecular orbitals of CHD and cZc-HT are frequently represented with the Mulliken symbols for irreducible representations of the C${}_2$ point group. Thus, the molecular orbitals 1${\sigma }^2$ 1${\pi }^2$ 2${\pi }^2$, 1${\sigma }^2$ 1${\pi }^2$ 2${\pi }^1$ 3${\pi }^1$ and 1${\sigma }^2$ 1${\pi }^2$ 3${\pi }^2$ of CHD are represented as 1${}^1$A, 1${}^1$B, and 2${}^1$A, respectively and in the same way, the molecular orbitals 1${\pi }^2$ 2${\pi }^2$ 3${\pi }^2$, 1${\pi }^2$ 2${\pi }^2$ 3${\pi }^1$ 4${\pi }^1$ and 1${\pi }^2$ 2${\pi }^2$ 4${\pi }^2$ of cZc-HT are represented as 1${}^1$A, 1${}^1$B, and 2${}^1$A, respectively. Now lets consider the following case, incoming photon generates population to the first excited state (1${}^1$B) of CHD, then one we can think of a population transfer to the first excited state of cZc-HT (1${}^1$B), as it was shown in the analysis of the individual orbital in Figure 2, but as mentioned above, this behaviour has not been observed. An alternative doorway can be opened if there is a crossing (i.e. internal conversion driven by the nonadiabatic interaction) between the bright 1${}^1$B and the dark 2${}^1$A states of CHD, through which the population of the bright state is transferred to the dark state. The inclusion of the dark 2${}^1$A state, proposed by van der Lugt and Oosterhoff [24], allows an internal conversion between the blue line and the dashed line (see Figure 3). Subsequently, a second internal conversion occurs from 2${}^1$A to 1${}^1$A, that results in the population transfer to the ground state of cZc-HT or CHD without the need to pass through the first excited state of cZc-HT, 1${}^1$B.

2.2. Spin-Flip TDDFT

The main application of TDDFT is the calculation of excited states which can be accessed through linear-response (LR) theory. In LR-TDDFT [43], the TD perturbation must be small enough such that the perturbation does not change completely the ground state of the system. Within LR-TDDFT, the non-Hermitian eigenvalue equation (1) $\left(\begin{array}{cc}A& B\\ B^{\ast }& A^{\ast }\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right)=\mathrm{\Omega }\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right)\left(\begin{array}{c}X\\ Y\end{array}\right)$(1) is solved for obtaining the excitation energies Ω and transition amplitudes X and Y. Here, the coupling matrices A and B are: (2) $A_{\mathrm{ia},\mathrm{jb}}=({ϵ}_a-{ϵ}_i){\delta }_{\mathrm{ab}}{\delta }_{\mathrm{ij}}+⟨\mathrm{ib}|\mathrm{aj}⟩-C_x⟨\mathrm{ib}|\mathrm{ja}⟩+⟨\mathrm{ib}|f^{\mathrm{xc}}|\mathrm{aj}⟩$(2) and (3) $B_{\mathrm{ia},\mathrm{jb}}=⟨\mathrm{ij}|\mathrm{ab}⟩-C_x⟨\mathrm{ij}|\mathrm{ba}⟩+⟨\mathrm{ij}|f^{\mathrm{xc}}|\mathrm{ab}⟩,$(3) in which (4) $⟨\mathrm{pq}|\mathrm{rs}⟩=\iint \mathrm{d}{r}_1\mathrm{d}{r}_2{\varphi }_p(r_1){\varphi }_q(r_2)\frac{1}{r_{12}}{\varphi }_r(r_1){\varphi }_s(r_2)$(4) and (5) $\begin{array}{rl}& ⟨\mathrm{pq}|f^{\mathrm{xc}}|\mathrm{rs}⟩\\ & =\iint \mathrm{d}{r}_1\mathrm{d}{r}_2{\varphi }_p(r_1){\varphi }_q(r_2)f^{\mathrm{xc}}(r_1,r_2){\varphi }_r(r_1){\varphi }_s(r_2)\end{array}$(5) where indexes $i,j\cdots$ and $a,b\cdots$ refer to occupied and unoccupied spin-orbitals, respectively, and $C_x$ is the fraction of Hartree-Fock (HF) exchange in the exchange-correlation (xc) functional $f^{\mathrm{xc}}$. The latter is the xc kernel which is given by the functional derivative of the xc energy. In SF-TDDFT method, that was introduced by Shao et al. [31] in 2003, a high-spin ($M_s=1$) triplet state, with two unpaired α-electron, is chosen as a reference state. Note that for conventional LR-TDDFT calculations, only $\mathrm{\alpha \alpha }$ and $\mathrm{\beta \beta }$ blocks are considered. In SF-TDDFT, the target states must be $M_s=0$, thus only $\mathrm{\alpha \beta }$ blocks are taken into account. Using a collinear exchange-correlation functional kernel, i.e. $⟨\mathrm{ij}|f^{\mathrm{xc}}|\mathrm{ab}⟩=0$, and considering the orthogonality between the occupied α and β spin-orbitals, the coupling matrices equations, Equation 2(2) $A_{\mathrm{ia},\mathrm{jb}}=({ϵ}_a-{ϵ}_i){\delta }_{\mathrm{ab}}{\delta }_{\mathrm{ij}}+⟨\mathrm{ib}|\mathrm{aj}⟩-C_x⟨\mathrm{ib}|\mathrm{ja}⟩+⟨\mathrm{ib}|f^{\mathrm{xc}}|\mathrm{aj}⟩$(2) and Equation 3(3) $B_{\mathrm{ia},\mathrm{jb}}=⟨\mathrm{ij}|\mathrm{ab}⟩-C_x⟨\mathrm{ij}|\mathrm{ba}⟩+⟨\mathrm{ij}|f^{\mathrm{xc}}|\mathrm{ab}⟩,$(3) , reduce to the following expressions (6) $A_{\mathrm{ia},\mathrm{jb}}=({ϵ}_a-{ϵ}_i){\delta }_{\mathrm{ab}}{\delta }_{\mathrm{ij}}-C_x⟨\mathrm{ib}|\mathrm{ja}⟩$(6) and (7) $B_{\mathrm{ia},\mathrm{jb}}=0.$(7) Note that Equation 7(7) $B_{\mathrm{ia},\mathrm{jb}}=0.$(7) is a result of using a collinear kernel which is similar to set B = 0 within the Tamm-Dancoff approximation (TDA)[44] in the LR-TDDFT. Thus, the SF-TDDFT eigenvalue equation to be solved is given by (8) $\mathrm{AX}=\mathrm{\Omega X}.$(8) SF-TDDFT is a very efficient method for describing excitation energies, topology around conical intersections, in particular those involving ground state, and excited-state reaction pathways. However, this method has one big drawback that is spin contamination and an effective state-tracking algorithm is required whenever SF-TDDFT is used. A straightforward way to do this, is to monitor the change in the excited-state transition density. An example, used in this work, is the state-tracking algorithm that is based on the overlap of the attachment/detachment densities at successive steps that can avoid any problems that may be introduced by sign changes in the orbitals or configuration-interaction coefficients. It has been shown by Herbert et al. [45] that although a state-tracking algorithm, in principle improve the performance of SF-TDDFT for CIs or minimum energy crossing points (MECP) optimisation, the spin adopted (SA)SF-TDDFT, being free of spin contamination and just as computationally efficient as SF-TDDFT, is a safer remedy.

2.3. Minimum-Energy crossing point optimisation

In nonadiabatic process, the crossing between two electronic states (I and J) along a $N_{\mathrm{int}}$-2-dimensional seam ($N_{\mathrm{int}}$ is the internal degrees of freedom) has a special significance [46,47] since the MECPs along this seam are one of the factors controlling the excited state nonadiabatic deactivation pathways [46,48]. The degeneracy between two states I and J is called conical intersection (CI) because the form of the two PESs around of seam has a conical shape, that arises because the degeneracy between the two states, I and J, is lifted along two orthogonal vectors; the difference gradient (DG) vector (9) $g^{\mathrm{IJ}}={\hat{\mathrm{\nabla }}}_R(E_I-E_J),$(9) which points in the direction of maximal energy splitting and the nonadiabatic coupling (NAC) vector (10) $h^{\mathrm{IJ}}=⟨{\mathrm{\Psi }}_I|{\hat{\mathrm{\nabla }}}_R|{\mathrm{\Psi }}_J⟩,$(10) which points in the direction of maximal nonadiabatic interaction (see Figure 4). These two vectors describe the 2-dimensional branching space defining the CI. For any electronic structure method that has analytic excited-state gradients, the $g^{\mathrm{IJ}}$ vector is already available, while computing $h^{\mathrm{IJ}}$ vector is not an easy task, since it requires first-derivative of the electronic wave function with respect to the nuclear coordinates. Although NAC is needed for optimisation of MECPs, there exist efficient MECPs optimisation algorithms that do not necessarily need calculation of NAC vector [49–51]. Recently, analytic derivative couplings for SF-TDDFT has been reported [52], thus here we rely on an optimisation algorithm that make use of both $g^{\mathrm{IJ}}$ and $h^{\mathrm{IJ}}$ vectors, i.e. the projected-gradient algorithm of Bearpark et al. [49] which is implemented in Q-CHEM [53]. This algorithm consists of minimising along the following gradient: (11) $G=2(E_J-E_I)x+P{G}_{\mathrm{mean}},$(11) where (12) $x=\frac{g^{\mathrm{IJ}}}{||g^{\mathrm{IJ}}||}$(12) is the normalised DG vector, (13) $G_{\mathrm{mean}}=\frac{1}{2}{\hat{\mathrm{\nabla }}}_R(E_I+E_J)$(13) is the mean energy gradient for states I and J, and (14) $P=1-x{x}^T-y{y}^T$(14) is a projector operator onto the seam space where (15) $y=\frac{(1-x{x}^T)h^{\mathrm{IJ}}}{||(1-x{x}^T)h^{\mathrm{IJ}}||}.$(15) It has been shown that the projected gradient algorithm, that uses both $g^{\mathrm{IJ}}$ and $h^{\mathrm{IJ}}$, reduces the number of iterations necessary to converge the MECP [52,54]. This fact makes the projected gradient algorithm more efficient than those algorithms that only use the vector $g^{\mathrm{IJ}}$.

Figure 4. Schematic representation of a CI between two PESs in the $N_{\mathrm{int}}$-2-dimensional space. Two orthogonal vectors define the 2-dimesional branching space: g which points in the direction of maximal energy splitting and h which points in the direction of maximal nonadiabatic interaction.

3. Computational details

Ground state geometry optimisation of CHD and cZc-HT were performed at ωB97X-D/cc-pVDZ level of theory [55,56], including Grimme's dispersion correction [57]. Subsequent frequency calculations were undertaken to confirm the nature of the minima. The vertical excitation energies of the first (S${}_1$) and the second (S${}_2$) excited state of CHD and cZc-HT and their optimised geometries were obtained using SF-TDDFT/cc-pVDZ in combination with the ωB97X-D, B5050LYP, and BHHLYP xc functionals, for bench-marking purpose, from which BHHLYP was finally selected. Relaxed potential energy surface (PES) scans of S${}_0$, S${}_1$ and S${}_2$ were carried out, along C${}_1$–C${}_6$ bond distance, at SF-TDDFT(BHHLYP)/cc-pVDZ level of theory. The geometries of all identified minima and transition state structure have been individually optimised. Furthermore, the nature of the transition state has been confirmed using the intrinsic reaction coordinate (IRC) algorithm, which is essentially a series of steepest descent steps going downhill from the transition state to the adjacent educt and product states. The MECPs between S${}_2$/S${}_1$ and S${}_1$/S${}_0$ were located using SF-TDDFT(BHHLYP)/cc-pVDZ. Finally, natural transitions orbitals (NTOs) of important geometries were computed using SF-TDDFT(BHHLYP)/cc-pVDZ. All the calculations were performed using Q-Chem 5.1.2 package [53]. Cartesian coordinates of all the relevant structures are given in the Supporting Information (SI).

4. Results and discussions

4.1. Stationary geometries

Ground state optimised geometries of CHD and cZc-HT molecules are illustrated in Figures 5 (a) and (b), respectively. The distance between C${}_1$ and C${}_6$ (reaction coordinate here) is 1.53 Å for the closed (CHD) and 3.44 Å for the open (cZc-HT) ring. Both geometries of CHD and cZc-HT show C${}_2$ symmetry in their S${}_0$ minimum. C${}_2$ has two irreducible representations, A (properties are symmetric to both the identity operation E and a rotation 180 degree rotation around C${}_2$ axis) and B (properties are symmetric with respect to the identity operation E). The ground state optimised geometries are in good agreement with previous works [21,22].

Figure 5. Optimised geometries with DFT(ωB97X-D)/cc-pVDZ: (a) ground state CHD, (b) ground state cZc-HT. Optimised geometries with SF-TDDFT(BHHLYP)/cc-pVDZ: (c) CI${}_{S_2/S_1}$ (from CHD), (d) CI${}_{S_1/S_0}$, (e) CI${}_{S_2/S_1}^{\mathrm{\prime }}$ (from cZc-HT), (f) S${}_{1min}$, (g) S${}_{2min}$ and (h) S${}_0$ transition state. Distances are in Å.

The S${}_1$ geometry optimisations of CHD and cZc-HT molecules both converged to the same local minimum with an energy of 3.92 eV, symmetry of C${}_1$, and the C${}_1$–C${}_6$ distance of being 2.16 Å (see Figure 5 (f)). Please note that the C${}_1$–C${}_6$ bond length is elongated (by 0.63 Å) and shortened (by 1.28 Å) compared with CHD-S${}_{0min}$ and cZc-HT-S${}_{0min}$, respectively. Such structural changes make it evident that vibrational relaxation on S${}_1$ leads to the ring opening of the CHD. The energy of S${}_{1min}$ is in a reasonable agreement with previously reported MS-CASPT2 results (4.40 eV), with the C${}_1$–C${}_6$ distance being larger, ≈ 0.5 Å  compared with the corresponding distance using MS-CASPT2 [18], and in both methodologies $S_1$ vibrational relaxation leads to the opened ring structure. Similarly, the S${}_2$ geometry optimisations of the CHD and cZc-HT molecules also converged to the same local minimum with an energy of 4.34 eV, symmetry of C${}_2$, and the C${}_1$–C${}_6$ distance of being 2.23 Å (see Figure 5 (g)). It must be noted that although vibrational relaxation on the S${}_2$ leads to the ring opening, but unlike S${}_{1min}$, the S${}_{2min}$ has C${}_2$ symmetry similar to CHD-S${}_{0min}$ and cZc-HT-S${}_{0min}$. To the best of our knowledge, no theoretical and experimental results have been reported for the S${}_{2min}$. It must be noted that unconstrained geometry optimisation of the S${}_1$ state leads directly to the crossing region between the S${}_0$ and S${}_1$ states, while unconstrained geometry optimisation of the S${}_2$ state leads to the crossing region between S${}_1$ and S${}_2$ states. Both of these points will be discussed later. Our results shows that interconversion between CHD and cZc-HT in the S${}_0$ proceed via a transition state with C${}_2$ symmetry, and the C${}_1$–C${}_6$ distance of being 2.0 Å (see Figure 5 (h)). The energy barrier for this process amounts to 3.4 eV, indicating that thermal electrocyclic reaction can be excluded from being relevant in the CHD to cZc-HT interconversion.

4.2. Vertical excitation energies

Table 1 summarises the vertical excitation energies (VEE) for the two lowest singlet excited states of the CHD and cZc-HT molecules calculated with TDDFT (various functionals) as well as experimental and previous theoretical studies using multi reference wave function-based methods. Recently, Polyak et al. [22] have performed a comprehensive theoretical study on the CHD photodissociation using the XMS-CASPT2 method. They reported 5.13 eV for the S${}_1$ VEE and 6.28 eV for the S${}_2$ VEE, that are very close to the experimental results [5] (see Table 1). Comparing the performance of SF-TDDFT with wave function-based methods listed in Table 1, it is apparent that functionals containing larger amount of HF exchange (>50%), perform in a closer agreement (difference of $\approx 0.1$ eV for S${}_1$ and $\approx 0.4$ eV for S${}_2$) with the corresponding value obtained at XMS-CASPT2/cc-pVDZ level of theory. For all level of theories, the order of states remains essentially the same; the S${}_1$ state with a 1${}^1$B symmetry being the bright $\pi {\pi }^{\ast }$ state and the S${}_2$ with the 2${}^1$A symmetry being the dark $\pi {\pi }^{\ast }$ state.

Table 1. Vertical excitation energies in eV for the two lowest singlet excited states of CHD and cZc-HT.

		CHD/States		cZc-HT/States
Level of theory		S${}_1$	S${}_2$	S${}_1$	S${}_2$
SF-TDDFT	ωB97X-D/cc-pVDZ	4.67	6.05	5.70	6.60
	BHHLYP/cc-pVDZ	5.26	5.89	6.23	6.69
	B5050LYP/cc-pVDZ	5.26	5.90	6.26	6.72
DFT & Mul. Ref. Methods	MS-CASPT2 [18]	5.27	6.39	–	–
	SA3-CASSCF [18]	7.44	6.55	–	–
	SSR-ωPBEh/6-31G* [21]	4.78	–	6.06	–
	XMS-CASPT2/cc-pVDZ [22]	5.13	6.28	–	–
Experiment [5]		4.94	6.30	–	–

Additionally, for the sake of completeness, we calculated the VEE of cZc-HT (see Table 1). Similar to the CHD molecule, all functionals shown that the S${}_1$ state with a 1${}^1$B symmetry is the bright $\pi {\pi }^{\ast }$ state and the S${}_2$ with the 2${}^1$A symmetry is the dark $\pi {\pi }^{\ast }$ state (The state-averaged NTO involved in the transitions are shown in Supporting Information Figure 1). To the best of our knowledge, no experimental results is available for the cZc-HT molecule. We observed almost no difference in the excitation energies obtained using BHHLYP and B5050LYP. Earlier studies [31–33] have revealed that BHHLYP is a suitable functional when describing excited sate and MECP optimisations within SF-TDDF methodology. Follwoing these conclusions, we will restrict ourselves to BHHLYP functional for further discussion.

4.3. Minimum energy crossing points

As well known [46–48], the excited state nonadiabatic dynamics is essentially controlled by the energies of the minima of the various surfaces, the various conical intersection seams, and MECPs along these seams. In this section, we discuss two MECPs that are known to play significant roles in actual photochemical interconversion of CHD to cZc-HT, namely, the CIs between S${}_2$ and S${}_1$ and S${}_1$ and S${}_0$, denoted as CI${}_{S_2/S_1}$ and CI${}_{S_1/S_0}$, respectively. Optimised structures of these geometries are shown in Figure 5 and their relatives energies to S${}_{0min}$ of CHD are listed in Table 2. Furthermore, the energies of the two states involved in the MECP optimisation, norm of derivative coupling (i.e. NACs) and the C${}_1$–C${}_6$ distance along the MECPs optimisation trajectories, are plotted in Figures 6 – 7.

Figure 6. (a) CI${}_{S2/S1}$: potential energies curves (i), norm of derivative coupling (ii) and distance between the reaction carbons C${}_1$–C${}_6$ (iii) along the S${}_2$/S${}_1$ MECP optimisation trajectory of the CHD. (b) CI${}_{S1/S0}$: potential energies curves (i), norm of derivative coupling (ii) and distance between the reaction carbons C${}_1$–C${}_6$ (iii) along the S${}_1$/S${}_0$ MECP optimisation trajectory of the CHD. The energies are relatives to the S${}_{0min}$ energy of the CHD.

Figure 7. (a) CI${}_{S2/S1}^{\mathrm{\prime }}$: potential energies curves (i), norm of derivative coupling (ii) and distance between the reaction carbons C${}_1$–C${}_6$ (iii) along the S${}_2$/S${}_1$ MECP optimisation trajectory of the cZc-HT. (b) CI${}_{S1/S0}$: potential energies curves (i), norm of derivative coupling (ii) and distance between the reaction carbons C${}_1$–C${}_6$ (iii) along the S${}_1$/S${}_0$ MECP optimisation trajectory of the cZc-HT. The energies are relatives to the S${}_{0min}$ energy of the CHD.

Table 2. C${}_1$–C${}_6$ bond distances and energies at different level of theories. The distances are in Å and after the slash, the energies relative to S${}_{0min}$ energy of CHD are in eV.

Method/basis set		C${}_1$ –C${}_6$/S${}_{1\mathrm{FC}}$	C${}_1$ –C${}_6$/S${}_{1min}$	C${}_1$–C${}_6$/S${}_{2min}$	C${}_1$–C${}_6$/CI${}_{S_2/S_1}$	C${}_1$–C${}_6$/CI${}_{S_1/S_0}$	C${}_1$–C${}_6$(cZc-HT)/S${}_{1\mathrm{FC}}$
SF-BHHLYP/cc-pVDZ		1.53/5.26	2.16/3.92	2.23/4.34	2.03/4.63	2.14/3.97	3.44/6.23
					2.36/4.42${}^a$
XMS-CASPT2/cc-pVDZ [22]		1.54/5.13	–	–	2.03/4.39	2.14/3.80	–
SSR-ωPBEh/6-31G* [21]		1.529/4.78	1.703/4.01	–	–	2.129/4.13	3.461/6.06
MS-CASPT2 [18]		1.533/5.27	1.727/4.40	–	1.986/4.51	2.137/3.64	–
SA3-CASSCF [18]		1.533/7.44	2.067/6.14	–	1.862/6.39	2.338/4.66	–
Experiment [3,5]		1.534/4.94	–	–	–	–	–

${}^a$ CI${}_{S_2/S_1}^{\mathrm{\prime }}$

CI${}_{S_2/S_1}$. The CHD and cZc-HT S${}_{0min}$ geometries served as the starting point in both calculations searching for the MECPs along the crossing seam between the S${}_2$ and S${}_1$ states. Starting from CHD, the optimisations converged to a minimum, CI${}_{S_2/S_1}$, with an energy of 4.63 eV, symmetry of C${}_2$, and the C${}_1$–C${}_6$ distance of being 2.03 Å (see Table 2 and Figure 5 (c)). Interestingly, the nature of this structure is very similar to the S${}_{2min}$ with respect to the symmetry, energy (difference of 0.29 eV higher) as well as the C${}_1$–C${}_6$ distance (difference of 0.2 Å). The latter will indicate that the CI${}_{S_2/S_1}$ is located in the close vicinity of the S${}_{2min}$. The geometry of this MECP is in good agreement with previous works [18,22]. The energy of the CI${}_{S_2/S_1}$ is 0.63 eV lower than the energy of the CHD molecule at the FC point, according to the SF-BHHLYP/cc-pVDZ method; this energy difference is 0.74 eV for XMS-CASPT2/cc-pVDZ and 0.76 eV for MS-CASPT2 calculations (see Table 2). Thus, it is natural to assume that, for all these three methods, the CHD will undergo vibrational relaxation downhill from the FC region and reach the CI${}_{S_2/S_1}$. Starting from cZc-HT, the optimisations converged to a minimum, CI${}_{S_2/S_1}^{\mathrm{\prime }}$, with an energy of 4.42 eV, symmetry of C${}_2$, and the C${}_1$–C${}_6$ distance of being 2.36 Å (see Table 2 and Figure 5 (e)). Interestingly, the nature of this structure is also very similar to the S${}_{2min}$ with respect to the symmetry, energy (difference of 0.08 eV) as well as the C${}_1$–C${}_6$ distance (difference of 0.13 Å). We observe that the CI${}_{S_2/S_1}^{\mathrm{\prime }}$ is much closer to the S${}_{2min}$. The energy of the CI${}_{S_2/S_1}^{\mathrm{\prime }}$ is 1.81 eV lower than the energy of the cZc-HT molecule at the FC point, according to the SF-BHHLYP/cc-pVDZ method (see Table 2). Note that the CI${}_{S_2/S_1}^{\mathrm{\prime }}$ is only 0.21 eV lower than the CI${}_{S_2/S_1}$, and the S${}_{2min}$ is located between these two CIs, being closer to the former (see Table 2). To the best of our knowledge, the CI${}_{S_2/S_1}^{\mathrm{\prime }}$ has not been previously reported. In sum, we propose that after photoexcitation to S${}_1$, both the CHD and cZc-HT undergo vibrational relaxation to S${}_{1min}$ after passing through low-lying MECPs between S${}_2$ and S${}_1$, i.e. CI${}_{S_2/S_1}$ and CI${}_{S_2/S_1}^{\mathrm{\prime }}$, respectively. This will be discussed later in the context of the PESs for the overall the CHD to cZc-HT photochemical interconversion (see Figure 8).

Figure 8. Schematic representation of the PESs of CHD/cZc-HT photochemical interconversion. The reaction coordinate is the C${}_1$–C${}_6$ bond distance. The 1${}^1$A (ground state) is in blue, the state 1${}^1$B is in green and the state 2${}^1$A is in red. The curves 1${}^1$B and 2${}^1$A are a pictorial description that connects the important geometries computed with SF-BHHLYP/cc-pVDZ. The curve 1${}^1$A is a PES scan trough the C${}_1$–C${}_6$ bond distance computed with SF-BHHLYP/cc-pVDZ. The cones in purple represent the MEPCs. The arrows depict how the photochemical interconversion process follow after an absorption of a photon generates the population of the first excited state (1${}^1$B) of CHD.

CI${}_{S_1/S_0}$. The CHD and cZc-HT S${}_{0min}$ geometries served as the starting point in both calculations searching for the MECPs along the crossing seam between the S${}_1$ and S${}_0$ states. Interestingly, the optimisations converged, independent of starting from CHD or cZc-HT, to the same local minimum with an energy of 3.97 eV, symmetry of C${}_1$, and C${}_1$–C${}_6$ distance of being 2.14 Å (see Table 2 and Figures 5 (d)). Interestingly, the nature of this structure is very similar (almost identical) to the S${}_{1min}$ with respect to the symmetry, energy (difference of 0.05 eV higher) as well as the C${}_1$–C${}_6$ distance (difference of 0.02 Å). The latter speaks in favour of the fact that the CI${}_{S_1/S_0}$ is located extremely close to the S${}_{1min}$. The geometry of this MECP is in good agreement with previous works [18,21,22]. The energy of the CI${}_{S_1/S_0}$ is 0.66 eV lower than the energy of the CI${}_{S_2/S_1}$ geometry, according to the SF-BHHLYP/cc-pVDZ method; the corresponding energy difference is energy is 0.59 eV for XMS-CASPT2/cc-pVDZ and 0.87 eV for MS-CASPT2 calculations (see Table 2). Thus, it is natural to assume that, for all these three methods, the CHD will go downhill from the CI${}_{S_2/S_1}$ towards CI${}_{S_1/S_0}$, and from there towards either to the S${}_{1min}$, that is slightly lower in energy (≈0.05 eV), or to the S${}_0$ surface. This fact suggests that the relaxation process from CI${}_{S_2/S_1}$ to CI${}_{S_1/S_0}$ is accelerated, supporting the hypothesis formulated by Garavelli et al. [7] about the ultrashort relaxation lifetime reported after S${}_2$/S${}_1$ crossing [6], i.e. there is a higher probability of hitting a CI during motion along the bottom of S${}_{1min}$. The similar scenario holds if one considers relaxation starting from CI${}_{S_2/S_1}^{\mathrm{\prime }}$.

4.4. Potential energy surface

In this section, the PESs for the lowest three singlet states, S${}_0$, S${}_1$ and S${}_2$, along C${}_1$–C${}_6$ coordinate, are investigated in more detail, shedding further light on the CHD/cZc-HT photochemical interconversion. We use the S${}_0$, S${}_1$, S${}_2$ to indicate adiabatic states and the $1^{1}A$, $1^{1}B$ and $2^{1}A$ to refer diabatic states. As it is depicted in Figure 8, the the S${}_0$ PES has two minima corresponding to i) the CHD with C${}_1$–C${}_6$ distance of 1.53 Å (closed ring) and ii) the cZc-HT with C${}_1$–C${}_6$ distance of 3.44 Å (open ring). The energy barrier for connecting these two minima, along C${}_1$–C${}_6$ coordinate, is 3.4 eV (78.41 kcal/mol), via the transition state structure at 2.0 Å. This large barrier practically excludes the thermal electrocyclic interconversion as a plausible reaction, and supports the photochemical electrocyclic reaction being the operative pathway.

After absorption of a photon, the CHD is vertically excited, from its S${}_{0min}$, to the first (bright) singlet excited state, S${}_1$(1${}^1$B), with an energy of 5.26 eV. From the FC region, the CHD undergoes a vibrational relaxation, during which it hits the CI${}_{S2/S1}$ (at 2.03 Å and energy of 4.63 eV), via which internal conversion occurs. This behaviour is typical when the CI is accessible from the FC region without significant energy barriers that is the case here. From this critical point, the reaction pathway may bifurcate in two directions, depending on which character the wave function takes after the system leaves the crossing region. Two scenarios can be defined; 1) if the diabatic character changes from 1${}^1$B to 2${}^1$A, then the system will evolve directly toward the CI${}_{S1/S0}$ (at 2.14 Å with an energy of 3.97 eV) that triggers an ultra-fast internal conversion process and provides a funnel of fast access to the ground state, on which the system can evolve either to the S${}_{0min}$ of CHD or the S${}_{0min}$ of cZc-HT. 2) if the diabatic character does not change, system evolves towards the S${}_{2min}$(1${}^1$B), from where the low-lying CI${}_{S_2/S_1}^{\mathrm{\prime }}$ (at 2.36 Å and energy of 4.42 eV) is accessible without significant energy barriers (only 0.08 eV higher), keeping in mind the excess vibrational energy after photoexcitation. After internal conversion through the CI${}_{S_2/S_1}^{\mathrm{\prime }}$, the system will undergo a vibrational relaxation through which it slides down again hitting the CI${}_{S1/S0}$, that acts as again as a doorway for an ultrafast internal conversion to the ground state, on which the system can again evolve either to the S${}_{0min}$ of CHD or the S${}_{0min}$ of cZc-HT. To the best of our knowledge, the excited state deactivation pathway via the CI${}_{S_2/S_1}^{\mathrm{\prime }}$, has not been reported so far.

5. Conclusions

In the present work, we have applied SF-TDDFT to reinvestigate the photochemical interconversion between CHD and cZc-HT. Our results reveal that the SF-TDDFT/cc-pVDZ in combination with the BHHLYP functional was able to describe successfully the critical geometries of the ground, S${}_0$, and the first two excited-state, S${}_1$ and S${}_2$, potential energy surfaces, such as, various minima, transition state, minimum-energy crossing points between S${}_2$/S${}_1$ and S${}_1$/S${}_0$, being in good agreements with the corresponding structures obtained by multireference wave function methods [16,18,22] and a variant of DFT [21]. Additionally we have identified an additional MECP between S${}_1$ and S${}_2$, namely CI${}_{S_2/S_1}^{\mathrm{\prime }}$, that opens an alternative deactivation channel for the excited CHD molecule, and subsequently fully described the photochemical interconversion between CHD and cZc-HT. This suggest that SF-TDDFT could be a good low-cost method to study complex molecules that contain the CHD chromophore as their central backbone, such as dithienylethene molecular photo switches relevant in molecular electronic devices [11,58].

Acknowledgments

This work is part of Innovational Research Incentives Scheme Vidi 2017 with project number 016.Vidi.189.044, which is financed by the Dutch Research Council (NWO).

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is part of Innovational Research Incentives Scheme Vidi 2017 with project number 016.Vidi.189.044, which is financed by the Dutch Research Council.