The deactivation of π-stacked cytosine molecules following excitation by ultrashort laser pulses was studied using semiclassical dynamics simulations. Another deactivation channel was found to compete with a previously reported path that led to dimerization. For both pathways, the initial excited state was found to form a charge-separated neutral exciton state, which forms an excimer state by charge transfer. When the interbase distance becomes less than 3 Å, charge recombination occurs due to strong intermolecular interaction, ultimately leading to an avoided crossing. Results indicate that the C2–N1–C6–C5 and dihedral angles play a significant role in the vibronic coupling between the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO). Vibrational energy distribution determines the fate of the excimer at the avoided crossing. Higher-amplitude vibration of C5 or C6 atoms leads to a nonadiabatic transition to the electronic ground state (a photophysical pathway); otherwise, a chemical reaction leading to the formation of cyclobutane type dimer occurs as found in earlier studies. The S1 and S0 potential energy surfaces calculated at TD-DFT level and the simulated trajectories were found to be consistent with CASPT2 results.

1. Introduction

Electronic excitation of DNA nucleobases induced by solar ultraviolet light has been extensively studied since the 1960s because it can lead to the formation of harmful photoproducts such as cyclobutane pyrimidine dimers (CPDs) between adjacent pyrimidine bases within the same DNA strand [13]. The formation of CPDs is considered as the main cause of cell death, mutagenesis, and development of skin cancers [4]. However, the quantum yield of photolesion formation is very low, suggesting a highly efficient nonradiative deactivation of excited DNA to the ground state consistent with the photostability of nucleobases [3].

Electronically excited single nucleobases decay to the ground state primarily by ultrafast internal conversion in less than 1 picosecond. However, in oligo- and polynucleotides the decay of excited nucleobases is more complicated, (showing multiexponential processes with different time constants [514]) because of interbase interactions, including base stacking and base pairing, and base-solvent interactions [1518]. Longer-time decay observed in single-stranded oligonucleotides, where base paring is absent, is attributed to an intrastrand excimer/exciplex state resulting from π-stacking of two adjacent nucleobases [13, 1924].

The formation of a delocalized excimer state slows down ultrafast nonradiative decay [1924]. The ab initio calculations at CASPT2/CASSCF level [25] for two stacked adenines of different arrangements favor a mechanism involving two possible decay pathways. In this mechanism, unstacked or poorly stacked pairs of bases relax to the ground state by an ultrafast internal conversion through the conical intersection (CI) between the lowest excited state and ground state of the monomer. On the other hand, an excimer state formed between intrastrand stacked bases exhibits a longer decay time to the ground state because of an energy barrier along the path toward the CI of the monomer. It has been suggested [23] that the excimer state has a charge-transfer character and the state lifetime is correlated with the energy required to transfer an electron from one base to its stacked neighbor.

The formation of excimers is important not only to understand the distinct photophysics of oligonucleotides and DNA but also to account for the intrinsic and distinct photoinduced reaction of cytosine and thymine, which form CPDs. It has been proposed [14] that the singlet excimer state is a precursor to the photodimerization of DNA bases and the excited state dimerization reaction occurs in competition with an internal conversion to the electronic ground state of nucleobases. A recent theoretical calculation [26, 27] for cytosine photodimer formation at the CASSCF/CASPT2 level suggests that the nonadiabatic photodimerization of two stacked cytosine molecules proceeds through a CI involving the lowest singlet excited state and the ground state. This conical intersection, , occurs below the energy of excimer 1(CC). In order to reach from 1(CC), the system has to surmount an energy barrier of about 0.2 eV. The formation of the stable excimer along the lowest singlet excited-state for stacked cytosine molecules may make the formation of cyclobutane cytosine dimer (CC) less efficient compared to the formation of cyclobutane thymine dimer (TT) from a stacked pair of thymine molecules, which do not form a stable excimer. For the latter, the favored mechanism involves a barrierless nonadiabatic mechanism for the photodimerization along a singlet excited-state through a , which is the funnel for ultrafast nonadiabatic decay leading the formation of cyclobutane thymine dimer (TT). The leads to two possible processes of ultrafast internal conversion to the ground state. One pathway involves essentially a chemical change (dimerization), the other (photophysical) leads to two separated monomer bases.

In this paper we present a semiclassical dynamics simulation study for the laser-induced excitation and deactivation of two π-stacked cytosine molecules. The simulation results provide detailed information on the dynamics from photon excitation to deactivation and are expected to be helpful in understanding this process.

2. Methodology

The SERID method is used to carry out dynamics simulations. This technique is described in detail elsewhere [28, 29]. In this methodology, the valence electrons are calculated by the time-dependent Schrödinger equation while both the radiation field and the motion of the nuclei are treated classically. The one-electron states are obtained for each time step by solving the time-dependent Schrödinger equation in a nonorthogonal basis: where is the overlap matrix of the atomic orbitals. The laser pulse is characterized by the vector potential , which is coupled to the Hamiltonian via the time-dependent Peierls substitution [30]. Consider where is the Hamiltonian matrix element for basis functions and on atoms at and , respectively, and is the charge of the electron.

In SERID, forces acting on nuclei or ions are computed by the Ehrenfest equation: where is effective nuclear-nuclear repulsive potential and is the expectation value of the time-dependent Heisenberg operator for the α coordinate of the nucleus labeled by (with ). Equation (3) is obtained by neglecting the second and higher order terms of the quantum fluctuations in the exact Ehrenfest theorem.

The time-dependent Schrödinger equation (1) is solved by using a unitary algorithm obtained from the equation for the time evolution operator [31]. Equation (3) is numerically integrated with the velocity Verlet algorithm. A time step of 0.05 femtoseconds is used for this study and energy conservation was then found to hold better than 1 part in 106 in a one picosecond simulation at 298 K.

The present “Ehrenfest” principle is complementary to other methods based on different approximations, such as the full multiple spawning model developed by the Martínez group [32]. The limitation of this method is that the simulation trajectory moves along a path dominated by averaging over all the terms in the Born-Oppenheimer expansion: rather than following the time evolution of a single potential energy surface, which is approximately decoupled from all the others. (Here, and represent the sets of nuclear and electronic coordinates, resp., and the are eigenstates of the electronic Hamiltonian at fixed .) The strengths of the present approach include the retention of all of the nuclear degrees of freedom and incorporation of both the excitation due to a laser pulse and the subsequent deexcitation at an avoided crossing near a CI.

3. Results and Discussions

In a recent paper [33], we described simulations that suggest the formation of cyclobutane cytosine dimer (CC) from a π-stacked pair of cytosine molecules irradiated by an ultrashort laser pulse (4.1 eV, fluence = 83.62 J·m−2, fwhm = 25 fs). Formation of an exciton intermediate upon excitation, followed by excimer formation, was involved in this dimerization channel. This paper describes another deactivation pathway found from additional simulations involving the same initial configuration. A typical trajectory (frequency of 4.1 eV, laser duration of 50 fs and fluence of 77.19 J·m2) is detailed below.

The variations with time of the HOMO-1, HOMO, LUMO, and LUMO+1 energies and the time-dependent population of these four orbitals are shown in Figures 1(a) and 1(b), respectively. Figure 1(a) shows that there is an abrupt fall in the LUMO energy soon after the application of the laser pulse. An intersection, which named avoided crossing, between the HOMO and LUMO levels is induced by coupling of orbitals, with an energy gap of 0.02 eV, is found at 976 fs. Figure 1(b) shows that by the end of the laser pulse (at 50 fs), about 1.4 electrons are excited from the HOMO to the LUMO, implying electronic excitation of one cytosine molecule. The coupling between the HOMO and LUMO, observed at 976 fs in Figure 1(a), suggests electronic transition from the LUMO to the HOMO. This deexcitation ultimately brings the molecules to the electronic ground state. It can be seen from Figure 1(a) that, shortly after the coupling, both the LUMO and HOMO levels move toward their initial values. After the electronic transition, these two energy levels slightly fluctuate about constants values which are essentially the same as their initial values. The excited state lifetime of 920 fs is about 3~4 times that which is for nonadiabatic deactivation of cytosine monomer [34] and is consistent with dimerization [33].

The variations with time of the lengths of the C5–C6, , C5 and C6 bonds are shown in Figure 2. Both C5–C6 and are double bonds at the beginning of simulation. Starting at 1.38 Å, which is the length of a typical C–C double bond in a conjugated aromatic system, the C5–C6 bond length rapidly elongates to about 1.5 Å after the laser pulse is applied while the bond length remains essentially constant. This suggests excitation of the molecule, but not the molecule. After 1000 fs, both C5–C6 and bond lengths shorten to their initial values and remain at this length until the end of the simulation. In Figure 2(b), both C5 and C6 distances are more or less constant up to 600 fs. The C5 and C6 distances sharply fall at around 700 fs, both reaching a minimum value of 2.0 Å at about 1000 fs. The minima suggest that the CC is not formed in this trajectory. The structures inserted in Figure 2(b) represent the average geometries at different stages of the trajectory. The timescale from 100 to 700 fs is defined as stage I. In this stage, weaker intermolecular interaction is observed because of the relatively large distance (>3 Å) between the C and molecules. Excitation distorts the molecule (bottom); note that the unexcited molecule () remains planar. With the shortening intermolecular distance, the molecule deforms as the excimer evolves to stage II (900–1000 fs). After 1000 fs (stage III), the two cytosine molecules move far apart due to deactivation.

Figure 2(a), stage I, shows previously unobserved variations in the C5–C6 and bond lengths. The C5–C6 bond length shortens from 1.53 Å to 1.48 Å; the bond length stretches from 1.38 Å to 1.42 Å. These variations correlate with intermolecular charge transfer, as seen in Figure 3. Before 100 fs (Figure 3(b)), the relatively flat curve for the net charge on molecule suggests no intermolecular charge transfer. Thereafter, about 1.4 electrons transfer from to (Figure 3(a)), after about 900 fs (when interbase distance is 3.1 Å), charge recombination drives electrons back and the two molecules return to electrically neutral states at the avoided crossing (976 fs). When electrons migrate to , the increasing electronic density between C5 and C6 results in a reduced C5–C6 bond length because of the increasing C5–C6 bond order. Conversely, decreasing electronic density in due to the weakening of π bond is manifested by an increase in the bond length.

The variations in the valence electrons of the C5 and C6 atoms, which are most actively involved in the excitation-deactivation process of pyridines, are also shown in Figure 3; Figure 3(d) zooms in on the 0 to 100 fs region of Figure 3(c). The valence electrons on the atoms are calculated by projection of single electronic wave functions to the orbitals of these atoms. Because of bond polarities due heteroatoms in the aromatic nucleus, valence electrons number of carbon atoms usually fluctuate about the 4 instead of being exactly equal to 4. For example, the valence electrons number of C5 and C6 are 4.31 and 3.86, respectively. Excitation by laser results in redistribution of π electrons. Consequently, the covalence number for the C5 atom drops from 4.31 to 3.9, while increasing to an average of about 4.5 for the C6 atom during the first 100 fs. Since there is minimal electron transfer between the two cytosine molecules during this time period (Figure 3(b)), this indicates a charge separation within the excited molecule. The charge separated neutral excited state with a lifetime of 100 fs is very similar to the “Frenkel exciton states” [14] formed by the coupling of 1 states (which are localized on single bases) of proximal nucleobases proposed by Gustavsson and Markovitsi and coworkers [1518]. According to some experimental and theoretical research [23, 3537], exciton formation tends to precede the formation of an excimer/exciplex. The shortening of C5–C6 and elongation of , as shown in Figure 2(a), suggest an exciton state evolving to a charge transfer state (i.e., excimer) after 100 fs.

The variation of dihedral angles, C2–N1–C6–C5 and , which describe the twisting of the C5–C6 and bonds, respectively, is compared to the HOMO-LUMO energy gap in Figure 4. The dihedral C2–N1–C6–C5 decreases from 0° to −27° within 120 fs, then fluctuates around an average value of −20°. It indicates that the molecule has been deformed due to excitation. Greater torsion results in a smaller HOMO-LUMO energy gap and a higher electron affinity for , leading to an electron transfer from to . The C2–N1–C6–C5 angle reaches a maximum of −35° at 986 fs and returns to its initial value after that time. On the other hand, the dihedral varies little until 900 fs, then sharply rises to 15° at 986 fs. The variation of the dihedral angle suggests that interbase interaction leads to deformation of after 900 fs. The strong torsion of gives rise to nonadiabatic transition, which leads to the electronic ground state.

The geometry at the avoided crossing is shown in Figure 5. The C2–N1–C6–C5 and dihedral angles reach their maximum values at 976 fs. This strong torsion causes a nonadiabatic transition to electronic ground state. For the molecular geometry taken at 976 fs, the C5 and C6 distances are 2.13 and 2.19 Å, respectively. These compare favorably with 2.27 and 2.17 Å obtained by CASSCF/CASPT2 calculation [26] at the CI between the lowest excited singlet state and the ground state. This indicates that the avoided crossing found at 976 fs is, indeed, in close proximity to the CI. The distances are also very similar to the dimerization channel that we previously reported. We believe that the avoided crossing may be identical for both pathways.

It is worth noting that the average interbase distance is 3 Å at 900 fs as the dihedral angle reaches −30° and electrons transfer back. We obtained the same result in our previous work [33]. This suggests that 3 Å is a minimum value for both deactivation pathways. At larger distances, the excited system evolves to an excimer characterized by charge transfer; at shorter distances, the strong intermolecular interaction may change the redox potential of bases and result in charge recombination, which leads to two neutral molecules. We have previously studied the effect of interbase distance (from 3.4 to 4.4 Å) using SERID [38]. The results suggest that 4.2 Å is the upper limit for the deactivation pathways associated with π stacking.

Due to the substantial defects of the mean field theory used in SERID, the electronic states are not explicit. The potential energy surfaces for the S1 and S0 states were calculated by using long-range corrected time-dependent density functional (LRC-TD-DFT) theory [3942] for comparison and the results are summarized in Figure 6. The Cam-B3LYP function, which was reported to describe charge transfer state well [4345], was used in this research. In this calculation, energies of serial single points taken along the simulation trajectory were calculated. A CI with energy value of 2.5 eV was observed at 990 fs in the TDDFT potential energy surface, which is very close to avoided crossing in our simulation. The non-optimized initial ground state geometry results in a “quasi vertical excited energy” of 3.2 eV, which is much lower than CASPT2 result of 4.41 eV [26]. However, the TD-DFT and CASPT2 calculations still yield comparable values for the energy gap between the Frank-Condon point and CI (0.7 eV and 0.9 eV, resp.). Based on [26], we believe that an optimized ground state will be located at −1.21 eV level in the current coordinate system. Therefore, the energy of the CI will be 3.71 eV, which agrees well with 3.5 eV from CASPT2. In [26], the CASPT2 calculation also yields a minimum in S1 that is lower than the CI by 0.2 eV, as same as our TD-DFT results. Our simulation indicates that from the minimum, the (CC) excimer has to overcome a barrier of 0.9 eV to reach the funnel. The stability of the excimer may decrease the effectiveness of internal conversion in the singlet manifold compared to the thymine excimer. This explains why stacked thymines have lifetime of (about 600 fs [38]) which is shorter than stacked cytosine; internal conversion of the thymine excimer goes through a barrierless pathway. The geometry of the CI is shown in Figure 5. The C5 and C6 are 2.13 and 2.19 Å respectively, compare to 2.27 and 2.17 Å respectively in results of CASPT2. The similar geometries suggest that this semiclassical dynamics simulation is credible.

In our previous study [33], π-stacked cytosines could form CPD by laser pulse irradiation. In the dimerization mechanism, an evolution of “ultrafast exciton” to “charge transfer excimer” and ultimately leading to formation of CPD was reported. It is hard to understand that the almost identical approach in this work follows a photophysical deactivation mechanism. In order to distinguish dimerization and photophysical deactivation mechanism of π-stacked cytosines, we counted multiple trajectories simulations results. Different initial stacking orientations, with average intermolecular distance of 3.4 Å and dihedral C5–C6 of 36°, were used in multiple trajectories simulations. Simulations gave rise to both dimerization and photophysical deactivation pathway. To compare the vibrational degree of freedom, seven dimerization pathways and seven excimer pathways were chosen and several geometrical parameters with respect to central C5–C6 bond, namely, the distance of C6–N1 and C6–H6 and dihedral angles of C2–N1–C6–C5 and C4–C5–C6–H6, are presented in Figure 7. The red solid lines represent dimerization channels and blue dashed lines represent the photophysical deactivations channels., respectively. Excited state lifetimes of the all fourteen trajectories are approximately 350 fs because of the relatively short intermolecular distance [38]. All of the four aforementioned degrees of freedom have larger vibrational amplitudes in the photophysical channel before nonadiabatic transition. It suggests that vibration decreases possibility of bonding between C5 or C6 atoms.

Obviously, through the CI in PES, both evolution forward to photoproducts and back to reactants are possible. Therefore, two competitive channels exist in π-stacked cytosines system. Figure 8 summarizes the changes in molecular geometry for representative trajectories from the competing pathways. Shortly after application of the laser pulse (<100 fs), a neutral dipolar coupling excited state with ultrashort lifetime, an exciton, is formed. After 120 fs, the twisting of dihedral angle C2–N1–C6–C5 increases the electron affinity of and leads to a transfer of 1.4 electrons from to . As a result, an charge transfer excimer is formed. When two bases approach to within a distance of about 3 Å, charge recombination occurs; the stacked cytosines revert to neutral. The system evolves to an avoided crossing due to vibrational coupling of the HOMO and LUMO induced by maximal deformation of and molecules. At avoided crossing, vibrational energy distribution determines the deactivation fate.

4. Conclusion

This paper presents simulation results for the photophysical deactivation of the π-stacked cytosine system following excitation subjected by a laser pulse with 4.1 eV photon energy. Only one cytosine is excited at the beginning. Excitation initially leads to the formation of a charge separated neutral exciton state with a lifetime of 100 fs. An “excimer state” results from charge transfer that occurs when the interbase distance is 3.8 Å in this simulation. Charge recombination occurs when the interbase distance shortens to about 3 Å, leading to neutral stacked bases that evolve into an avoided crossing prior to deactivation. Geometries taken from avoided crossing are very similar to the S1/S0 CI obtained by CASSCT/CASPT2 [40, 41]. Torsional vibrations (C2–N1–C6–C5 and ) play a significant role in vibronic coupling between the HOMO and LUMO and lead to nonadiabatic transition of the molecules to the ground state.

The geometries of avoided crossing are almost the same as those for dimerization channel reported previously. It indicates there are two competitive de-exciton channels at avoided crossing. Higher amplitude vibrations of C5 or C6 atoms favor the photophysical deactivation pathway.

To compare with other theoretical results, S1 and S0 PESs were calculated at TD-DFT level along with simulation trajectory. The CI was found to be very near the avoided crossing obtained by SERID method. The TD-DFT calculation yields a 0.2 eV barrier between the minimum of S1 and CI, indicating a relatively stable excimer in π-stacked cytosines system, which lessens the effectiveness of internal conversion compared to thymine excimer, which has a barrierless pathway from S1 to CI.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This work is supported by the Joint Funds of the National Natural Science Foundation of China (Grant No. U1330138) and the Natural Science Foundation Project of CQ CSTC (cstc2011jjA00009) and Project of the Science Technology Foundation of Chongqing Education Committee, China (no. KJ120516 and KJ120521).