Exploring the Excited-State Nonadiabatic Effects in the Semisaturated Planar Tetracoordinated Carbon Molecule C₇H₄

Abstract

We theoretically study the nonadiabatic relaxation dynamics of low-lying singlet excited-states of semisaturated planar tetracoordinated carbon molecule, C${}_7$H${}_4$. This molecule possesses a stable C${}_{2\mathrm{v}}$ ground-state equilibrium geometry. The three low-lying singlet states, S${}_1$, S${}_2$ and S${}_3$, lie in the energy gap of about 1.2 eV. The potential energy surfaces constructed within the quadratic vibronic coupling formalism reveal multiple conical intersections in the Franck-Condon region. Upon photoexcitation to S${}_3$, the wavepacket decays rapidly to lower states via these conical intersections. We also observe the wavepacket transfer to S${}_3$ during the initial wavepacket evolution on lower states, suggesting the nonadiabatic behavior of photoexcited planar C${}_7$H${}_4$.

1. Introduction

The ongoing research in planar tetracoordinated carbon (ptC) for over five decades has opened a new era in the chemistry of carbon. Molecules containing ptC or even higher coordination have been successfully suggested and synthesized. The deviation from the long-established concept of tetrahedral tetracoordinate carbon by van’t Hoff [1] and Le Bel [2] has raised the curiosity of researchers. Designing systems containing a planar tetracoordinated carbon center has remained a challenging task. The lone pair of electrons on the central ptC and the electron-deficient 3-centered 2-electron bonds make ptC unstable. The concept of ptC was first introduced by Monkhorst in 1968 [3]. Subsequently, using a planar methane model, Hoffmann and coworkers proposed strategies for stabilizing planar tetracoordinated carbon arrangements by electronic effects as well as steric effects [4,5].

Using those strategies, a number of molecules were proposed theoretically [6,7,8,9,10,11,12,13,14,15,16], though only a few were identified in the laboratory [17,18,19,20,21]. Merino and coworkers have reported a series of semisaturated ptC candidates containing cyclic hydrocarbons, which were created by combining C${}_5^{2-}$ moiety with saturated hydrocarbon fragments [22]. These were the first semisaturated cycles, containing a ptC stabilized only by electronic factors. The authors analyzed various ground-state properties, and they inferred that the multicentric nature of the bonding within the C${}_5^{2-}$ skeleton and the resulting electron delocalization provides stability to these molecules.

Although several studies have focussed on the various ground-state properties of ptCs, the excited-state properties of these molecules remain unexplored. In the present work, we focus on elucidating the excited-state dynamics of the semisaturated ptC, C${}_7$H${}_4$ (tricyclo[4.1.0.0${}^{1,3}$]hept-2,6-diene-2,7-diyl) (cf., Figure 1), using combined electronic structure computations and quantum dynamics simulations. We generate the potential energy surfaces (PESs) of low-lying singlet excited electronic states within the quadratic vibronic coupling (QVC) framework. Subsequently, we perform quantum-mechanical wavepacket simulations within the well-established multidimensional configurational time-dependent Hartree (MCTDH) method. Finally, we analyze the electronic populations and reduced nuclear densities and various stationary points of PESs to gain insights into the excited-state relaxation decay channels of the ptC molecule, C${}_7$H${}_4$.

2. Theoretical Framework

2.1. Electronic Structure Calculations

We perform the optimization and frequency computations associated with the ground-state equilibrium geometry of C${}_7$H${}_4$ within the density functional theory (DFT) using gas-phase conditions. These computations use the B3LYP [23]/6-311++G(d,p) level of theory. For excited-state computations, we rely on the time-dependent variant of DFT (TDDFT). We obtain vertical excitation energies and corresponding oscillator strengths using the long-range corrected $\omega$B97XD functional (0.22 HF at short-range and 1.00 HF at long-range) [24] in combination with 6-311++G(d,p) basis set. All the calculations are done using Gaussian 16 program package [25].

To verify the suitability of the selected functional, we have done benchmarking calculations of excited-states with various quantum chemical methods. Accordingly, we employ time-dependent variants of B3LYP (containing 0.20 Hartree–Fock (HF) exchange), CAM-B3LYP [26] (comprises of 0.19 HF exchange at short-range, and 0.65 HF at long-range), LC-$\omega$PBE [27] (no HF exchange at short-range and 1.00 HF at long-range) and M06-2X (comprises of 0.54 HF exchange) [28]. We use the 6-311++G(d,p) basis set for these computations. Further, we also evaluate excited-state energies using post-HF methods: equation of motion coupled cluster with single and double excitations (EOM-CCSD) method [29,30], resolution-of-the-identity second-order approximate coupled-cluster singles and doubles (RI-CC2) method [31], and algebraic diagrammatic construction method to second-order (ADC(2)) method [32]. These computations use the correlation consistent polarized valence double zeta (cc-pVDZ) basis set. The geometry optimized at the B3LYP/6-311++G(d,p) level of theory was used for all the above computations. We employ Gaussian 16 software package for TDDFT and EOM-CCSD computations and TURBOMOLE 7.4 software [33] for ADC(2) and RI-CC2 computations.

2.2. Vibronic Hamiltonian

Based on the computed vertical excitation energies, we are interested in studying the excited-state dynamics of the first three singlet excited-states of C${}_7$H${}_4$. For this, we construct a 3 × 3 vibronic Hamiltonian based on the well-established QVC approach [34]. The Hamiltonian can be expressed as:

$$\begin{array}{cc}\hfill \mathcal{H}& =({\mathrm{T}}_N+{\mathrm{V}}_0){\mathbf{1}}_{\mathbf{3}}+\hfill \\ \hfill & \left(\begin{array}{ccc}{\mathrm{E}}_{{\mathrm{S}}_1}^0+{\displaystyle \sum _{i\in a_1}}{\kappa }_i^{\left({\mathrm{S}}_1\right)}{\mathrm{Q}}_i+\frac{1}{2}{\displaystyle \sum _{i\in a_1}}{\gamma }_i^{\left({\mathrm{S}}_1\right)}{{\mathrm{Q}}_i}^2& {\displaystyle \sum _{i\in a_1}}{\lambda }_i^{({\mathrm{S}}_1-{\mathrm{S}}_2)}{\mathrm{Q}}_i& {\displaystyle \sum _{j\in a_2}}{\lambda }_j^{({\mathrm{S}}_1-{\mathrm{S}}_3)}{\mathrm{Q}}_j\\ & {\mathrm{E}}_{{\mathrm{S}}_2}^0+{\displaystyle \sum _{i\in a_1}}{\kappa }_i^{\left({\mathrm{S}}_2\right)}{\mathrm{Q}}_i+\frac{1}{2}{\displaystyle \sum _{i\in a_1}}{\gamma }_i^{\left({\mathrm{S}}_2\right)}{{\mathrm{Q}}_i}^2& {\displaystyle \sum _{j\in a_2}}{\lambda }_j^{({\mathrm{S}}_2-{\mathrm{S}}_3)}{\mathrm{Q}}_j\\ h.c.& & {\mathrm{E}}_{{\mathrm{S}}_3}^0+{\displaystyle \sum _{i\in a_1}}{\kappa }_i^{\left({\mathrm{S}}_3\right)}{\mathrm{Q}}_i+\frac{1}{2}{\displaystyle \sum _{i\in a_1}}{\gamma }_i^{\left({\mathrm{S}}_3\right)}{{\mathrm{Q}}_i}^2\end{array}\right)\hfill \end{array}$$

(1)

Here, T${}_N$ and V${}_0$ are the ground-state kinetic energy and potential energy operators and are expressed in terms of the dimensionless normal coordinate (Q) within the harmonic approximation. E${}_{\mathrm{m}}^0$ represents the vertical excitation energy of the excited-states (for S${}_{\mathrm{m}}$, m = 1, 2 and 3). ${\kappa }_i$ and ${\gamma }_i$ denotes the linear and quadratic intrastate coupling parameters, respectively, along the totally symmetric vibrational mode, i. ${\lambda }_j$ is the interstate coupling parameter along the non-totally symmetric vibrational mode, j. The interstate coupling constants between states of same symmetry (S${}_1$ and S${}_2$) are computed along the totally symmetric vibrational mode, a${}_1$, and is represented as ${\lambda }_i$. These coupling parameters are evaluated using the following expressions:

$${\kappa }^{\left({\mathrm{S}}_{\mathrm{m}}\right)}=\frac{\partial V_{{\mathrm{S}}_{\mathrm{m}}}}{\partial {\mathrm{Q}}_i}{|}_{{\mathrm{Q}}_0},i\in a_1$$

(2)

$${\gamma }^{\left({\mathrm{S}}_{\mathrm{m}}\right)}=\frac{{\partial }^2{V}_{{\mathrm{S}}_{\mathrm{m}}}}{\partial {{\mathrm{Q}}_i}^2}{|}_{{\mathrm{Q}}_0},i\in a_1$$

(3)

$${\lambda }_{i/j}^{({\mathrm{S}}_{\mathrm{m}}-{\mathrm{S}}_{\mathrm{n}})}={\left[\frac{1}{8}\frac{{\partial }^2}{\partial {\mathrm{Q}}_{i/j}^2}{|V_{{\mathrm{S}}_{\mathrm{m}}}\left(\mathrm{Q}\right)-V_{{\mathrm{S}}_{\mathrm{n}}}\left(\mathrm{Q}\right)|}^2\right]}_{{\mathrm{Q}}_0}^{1/2},i\in a_1;j\in a_2$$

(4)

where, $V_{{\mathrm{S}}_{\mathrm{m}}}$ and $V_{{\mathrm{S}}_{\mathrm{n}}}$ are the adiabatic potential energies of two different electronic states, ${\mathrm{S}}_{\mathrm{n}}$ and ${\mathrm{S}}_{\mathrm{m}}$, respectively; Q${}_0$ represents the ground-state equilibrium geometry of the molecule.

2.3. Dynamics Simulations

We simulate the quantum nuclear dynamics of the singlet excited-states of C${}_7$H${}_4$ by using multiconfiguration time-dependent Hartree (MCTDH) approach [35,36,37] which is designed to solve the time-dependent Schrodinger equation for multidimensional dynamical systems. In this method, the nuclear wavefunction, $\mathsf{\Psi }$, of a system with f degrees of freedom (DOF) is expressed as:

$$\mathsf{\Psi }({\mathrm{Q}}_1,\dots ,{\mathrm{Q}}_f,t)=\sum _{i_1=1}^{n_1},\dots ,\sum _{i_f=1}^{n_f}{R}_{i_1,\dots ,i_f}\left(t\right)\prod _{k=1}^f{\varphi }_{ik}^{\left(k\right)}({\mathrm{Q}}_k,t).$$

(5)

where Q${}_1$, …, Q${}_f$ are the nuclear coordinates of vibrational modes. $R_{i_1,\dots ,i_f}$ and ${\varphi }_{ik}^{\left(k\right)}$ denote the MCTDH expansion coefficients and single-particle functions (SPFs), respectively. $n_k$ represents the number of SPFs to describe the k-th DOF.

To reduce the memory requirement of the above treatment, we adopt a “mode combination” technique in which SPFs that can describe a set of DOF are used. The nuclear wavefunction can be then rewritten as a multi-configuration over p generalized particles:

$$\mathsf{\Psi }({\mathrm{Q}}_1^{\prime },{\mathrm{Q}}_2^{\prime },\dots ,{\mathrm{Q}}_f^{\prime },t)=\sum _{j_1=1}^{n_1},\dots ,\sum _{j_p=1}^{n_p}{R}_{j_1,\dots ,j_p}\left(t\right){\varphi }_{j_1}^{\left(1\right)}({\mathrm{Q}}_1^{\prime },t),\dots ,{\varphi }_{j_p}^{\left(p\right)}({\mathrm{Q}}_p^{\prime },t)$$

(6)

with

$${\varphi }_j^{\left(k\right)}({\mathrm{Q}}_k^{\prime },t)={\varphi }_j^{\left(k\right)}({\mathrm{Q}}_1,{\mathrm{Q}}_2,\dots ,{\mathrm{Q}}_w,t)$$

(7)

where Q${}_k^{\prime }$ = (Q${}_1$, Q${}_2$, …, Q${}_w$) represents the multidimensional coordinate for mode k.

Out of the 28 degrees of freedom (27 vibrational modes and a set of electronic states) of C${}_7$H${}_4$, a total of 14 (9 a${}_1$ and 5 a${}_2$) modes were selected based on the excitation strength of the modes. It is to be noted that S${}_1$–S${}_3$ and S${}_2$–S${}_3$ couplings are along a${}_2$ modes and coupling between S${}_1$ and S${}_2$, belonging to the same symmetry, are along a${}_1$ modes. The nuclear wavepacket generated on S${}_3$ PES has been propagated for 300 fs with a time step of 1 fs. The diabatic electronic populations and nuclear densities are then extracted to study the internal conversion dynamics. The Heidelberg MCTDH code version 8.5 Revision 11 [38] is employed for these calculations. Details of MCTDH dynamics such as mode combination, primitive basis, and SPFs are given in Supplementary Materials.

3. Results and Discussion

3.1. PESs and Conical Intersections

We collect the vertical excitation energies, oscillator strengths and symmetries of low-lying singlet excited-states of C${}_7$H${}_4$ computed at different levels of theory in

Table 1. Vertical excitation energies of C${}_7$H${}_4$ computed using different levels of theory. 6-311++G(d,p) is the basis set used for the TDDFT methods and cc-pVDZ is used for the wavefunction based methods.

Methods	S${}_1$	S${}_2$	S${}_3$
TD-$\omega$B97XD	4.3417 (B${}_1$, 0.001)	4.8192 (B${}_1$, 0.043)	5.5606 (B${}_2$, 0.049)
TD-B3LYP	4.0221 (B${}_1$, 0.005)	4.5325 (B${}_1$, 0.032)	5.0136 (A${}_1$, 0.010)
TD-CAMB3LYP	4.3389 (B${}_1$, 0.001)	4.7960 (B${}_1$, 0.045)	5.5310 (B${}_2$, 0.045)
TD-LC-$\omega$PBE	4.5102 (B${}_1$, 0.000)	4.9985 (B${}_1$, 0.050)	5.7237 (B${}_2$, 0.037)
TD-M06-2X	4.2035 (B${}_1$, 0.000)	4.5731 (B${}_1$, 0.039)	5.4983 (A${}_1$, 0.004)
EOM-CCSD	4.3121 (B${}_1$, 0.002)	4.9658 (B${}_1$, 0.032)	5.5197 (B${}_2$, 0.031)
ADC(2)	4.3239 (B${}_1$, 0.003)	4.9202 (B${}_1$, 0.030)	5.4091 (B${}_2$ ,0.042)
RICC2	4.2634 (B${}_1$, 0.004)	4.8697 (B${}_1$, 0.027)	5.3844 (B${}_2$, 0.046)

. All computational methods yield the same electronic symmetries for S${}_1$ and S${}_2$. We also note that all methods predict B${}_2$ symmetry for S${}_3$ except TD-B3LYP and TD-M06-2X. The latter two methods show A${}_1$ symmetry for this state. Concerning oscillator strengths, all methods predict a smaller value for S${}_1$ than higher electronic states. Here, the important feature of Table 1 is the S${}_1$–S${}_3$ energy gap; all methods predict a gap of about 1.2 eV. As these states lie within this energy gap, one would expect the molecule to follow nonadiabatic behavior upon photoexcitation. As the vertical energies (as well as the oscillator strengths) of the excited-states computed at the $\omega$B97XD/6-311++G(d,p) level of theory match well with the computationally expensive wavefunction methods, we perform all calculations using this method.

To explore the nuclear dependence of these electronic states, we plot adiabatic potential energy profiles (solid lines) against in-plane ring deformation vibrational coordinates (Q${}_{15}$) in Figure 2. The ab initio energies (plus harmonic potential) are shown as filled circles in this figure. We observe that the potential energy profiles generated within the QVC approach reproduce the ab initio energies quite well, confirming the suitability of the latter approach to study the dynamics happening within the Franck-Condon (FC) region. We observe multiple curve crossings between the excited- states of interest. Such crossings also occur in the other vibrational coordinate space, however those data are not shown here for brevity. It should be mentioned that these crossings would form the seam of conical intersection in the multidimensional space. Figure 3 depicts the conical intersections associated with the adiabatic PESs of S${}_1$, S${}_2$, and S${}_3$ states along the coordinates of ring deformation (Q${}_{15}$) and C=C stretch (Q${}_{22}$) vibrations. These vibrational modes are represented in Figure 4.

We also evaluate two other important quantities: minimum energy conical intersection and equilibrium minimum of excited-states. These quantities are crucial to explain the wavepacket dynamics happening on the coupled S${}_1$–S${}_2$–S${}_3$ PESs. Following expressions are employed to evaluate those quantities: [39]

$$\mathrm{MECI}=\mathrm{a}+\frac{{(\mathrm{b}-\mathrm{c})}^2}{2\mathrm{d}}-\frac{1}{2}\sum _{i=1}^N\frac{{\sigma }_i^2}{{\omega }_i}$$

(8)

where,

$$\begin{array}{c}\mathrm{a}=\frac{{\mathrm{E}}_1^0+{\mathrm{E}}_2^0}{2}\end{array}$$

(9)

$$\begin{array}{c}\mathrm{c}=\frac{{\mathrm{E}}_2^0-{\mathrm{E}}_1^0}{2}\end{array}$$

(10)

$$\begin{array}{c}\mathrm{b}=\sum _{i=1}^N\frac{{\delta }_i{\sigma }_i}{{\omega }_i}\end{array}$$

(11)

$$\begin{array}{c}\mathrm{d}=\sum _{i=1}^N\frac{{\delta }_i^2}{{\omega }_i}\end{array}$$

(12)

$$\begin{array}{c}{\sigma }_i=\frac{{\kappa }_i^{\left(2\right)}+{\kappa }_i^{\left(1\right)}}{2}\end{array}$$

(13)

$$\begin{array}{c}{\delta }_i=\frac{{\kappa }_i^{\left(2\right)}-{\kappa }_i^{\left(1\right)}}{2}\end{array}$$

(14)

$${\mathrm{E}}_m^{\min }={\mathrm{E}}_m^0-\frac{1}{2}\sum _{i=1}^N\frac{{\kappa }_i^2}{{\omega }_i}$$

(15)

where ${\omega }_i$ is the harmonic frequency of the ith totally symmetric vibrational mode and N is the total number of such modes.

The molecule might display unexpected nonradiative decay dynamics due to the presence of accessible conical intersections within the FC region (cf.,

Table 2. Relevant stationary points (state minima and MECI) on the coupled S${}_1$–S${}_2$–S${}_3$ PESs of C${}_7$H${}_4$.

Stationary Point	Energy (eV)
S${}_1^{\min }$	3.6039
S${}_2^{\min }$	4.7940
S${}_3^{\min }$	4.3006
S${}_1$/S${}_2$	4.7970
S${}_1$/S${}_3$	5.6376
S${}_2$/S${}_3$	5.0552

). For instance, the S${}_2$–S${}_3$ conical intersection lies ∼0.5 eV below the FC point of S${}_3$ (cf., Table 1 and Table 2). Hence, upon excitation to S${}_3$, the molecule can, in principle, decay rapidly via this conical intersection. We expect this intersection point and the S${}_1$–S${}_3$ conical intersection, that lies slightly above (∼0.1 eV) the FC point of S${}_3$, might play an important role in the decay dynamics.

3.2. Singlet Dynamics

To investigate the fate of C${}_7$H${}_4$ in the higher excited singlet states, we perform dynamics simulations by launching the initial wavepacket at the FC point of the “bright” S${}_3$. Figure 5a collects the time-dependent electronic population profiles obtained from this wavepacket propagation calculation. We observe a rapid population transfer from S${}_3$ to S${}_2$. The wavepacket evolving on S${}_3$ could easily access the S${}_2$/S${}_3$ MECI (∼5.06 eV, cf., Table 2) and thus promote rapid nonadiabatic population transfer.

We note that ∼70% depopulation of S${}_3$ happens within the first 40 fs of propagation time. We observe a sharp rise in the S${}_1$ population during this period. This sharp rise would emerge from the wavepacket decay via S${}_1$–S${}_3$ and S${}_1$–S${}_2$ conical intersections. Interestingly, the populations of the involved states remain almost unchanged after 50 fs. These observations suggest the remnant wavepacket localization at the respective equilibrium minimum of involved electronic states after the early nonadiabatic wavepacket transfer.

To further validate the nonadiabatic population transfer, we plot the reduced nuclear densities associated with the ring deformation mode, i.e., Q${}_{15}$, for three electronic states in Figure 6. The density would be maximum on S${}_3$ at t=0 fs as the wavepacket propagation starts on this state. We observe the rapid density reduction within 20 fs, and after that, the density remains unchanged until the end of propagation time. For S${}_2$, the density accumulates rapidly after a few femtoseconds and reaches a maximum at about 20 fs. After that, the density reduces up to 50 fs and remains unchanged till 300 fs. For S${}_1$, we observe a gradual increase of nuclear density, reaching a maximum value at about 20 fs. After that, the nuclear density on S${}_1$ would remain unchanged up to 300 fs. From these observations, we infer that the wavepacket decay to S${}_1$ would occur within a few tens of femtoseconds after photoexcitation to S${}_3$.

Next, we launch the initial wavepacket at the FC point of S${}_2$ and S${}_1$ to explore the decay dynamics of individual electronic states. The accessible S${}_1$/S${}_2$ MECI (∼4.80 eV, cf., Table 2) combined with the high interstate coupling via a${}_1$ modes results in an extremely rapid population transfer from S${}_2$ to S${}_1$ (cf., Figure 5b). The S${}_1$ population reaches an asymptotic value of ∼0.6 at the end of 300 fs. A minor population transfer to S${}_3$ can also be seen, attributed to the strong nonadiabatic coupling between S${}_2$ and S${}_3$ states.

The wavepacket propagated on the S${}_1$ PES shows an appreciable decay (∼40%) to S${}_2$ till the initial 20 fs and further, remains constant (cf., Figure 5c). Such population profile suggests that the molecule might be trapped at the equilibrium minimum of S${}_1$, which is ∼0.7 eV lower in energy from the S${}_1$ FC point. The molecule could no longer access any conical intersections, resulting in no noticeable population transfer. Similar profile is also observed in the variation of nuclear densities (along Q${}_{15}$ vibrational mode) across the excited-states with the wavepacket evolving on S${}_1$ (cf., Figure S1 in Supplementary Materials). It is noted that no state shows 100% decay within 300 fs. So dynamics simulations need to be conducted for a longer duration for the complete decay of the excited-states.

We note that the QVC model employed in this study would be helpful to gain insights into the FC dynamics. However, concerning the fluxional behavior of molecules containing ptC, it is necessary to adopt a theoretical model that involves possible excited-state dissociation pathways. Identifying such dissociative pathways through experimental investigations can be complicated due to nonadiabatic events occurring in the FC region. However, one could rely on the high-level on-the-fly dynamics simulations to study the photoproducts of ptC molecules.

4. Conclusions

In the present work, we provide vital insights on the nonadiabatic relaxation dynamics of the semisaturated ptC molecule, C${}_7$H${}_4$. PESs constructed by employing QVC formalism revealed multiple conical intersections among the three singlet excited-states, S${}_1$, S${}_2$, and S${}_3$. The molecule excited to the optically “bright” S${}_3$ returns to S${}_1$ on an ultrafast timescale via these conical intersections. Similarly, these crossing points also promote the wavepacket transfer to higher states (S${}_2$ and S${}_3$) while the wavepacket is initially on S${}_1$, demonstrating the nonadiabatic behavior of C${}_7$H${}_4$ irrespective of excitation energy. Although the present study has provided essential features of excited-state decay dynamics, a detailed quantum dynamics study using the multi-reference configuration interaction (MRCI) or multiconfigurational quasidegenerate perturbation theory (MCQDPT) PESs is planned to provide an accurate picture of the relaxation pathways of C${}_7$H${}_4$.