Article Cite This: J. Chem. Theory Comput. 2019, 15, 3730−3742
pubs.acs.org/JCTC
Curious Case of 2‑Selenouracil: Efficient Population of Triplet States and Yet Photostable Sebastian Mai,* Anna-Patricia Wolf, and Leticia Gonzaĺ ez* Institute of Theoretical Chemistry, Faculty of Chemistry, University of Vienna, Währinger Straße 17, 1090 Vienna, Austria
Downloaded via UNIV GRONINGEN on August 7, 2019 at 12:40:00 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.
S Supporting Information *
ABSTRACT: Excited-state MS-CASPT2 and ADC(2) quantum chemical calculations and nonadiabatic dynamics simulations show that 2-selenouracil is able to both efficiently populate and depopulate reactive triplet states in an ultrashort time scale. Thus, the heavier homologue of 2-thiouracil unites the ultrafast, high-yield intersystem crossing of 2-thiouracil with the short excitedstate lifetime and photostability of the parent nucleobase uraciltwo properties that have been traditionally thought to be diametrically opposed. Remarkably, while the S2 → S1 → T2 → T1 deactivation dynamics of 2-selenouracil is analogous to that of 2-thiouracil, the calculations show that the triplet lifetime of 2-selenouracil should decrease by up to 3 orders of magnitude in comparison to that 2-thiouracil, possibly down to the few-picosecond time scale. The main reasons for this decrease are the lack of a second T1 minimum, the enhanced spin−orbit coupling, and the reduction of the energy barrier to access the T1/S0 crossingin particular in aqueous solutioncompared to 2-thiouracil. Such unusual photophysical properties, together with its significant red-shifted absorption spectrum, could make 2-selenouracil a useful specialized photosensitizer for photodynamical therapy. ground state,23,24 a mechanism that makes these bases photostable. On the contrary, thio-nucleobases quickly undergo intersystem crossing (ISC) that produces reactive, diradicaloid triplet states. These triplet states can cause severe damage to biological tissue,25 which can either be exploited in photodynamic therapy26 or might become an unwanted side effect. Hence, understanding the photoreactivity of nucleobase analogues is critical to design their biological and medicinal applications. Following this motivation, the past few years have seen a rise in research about photoreactivity of nucleobase analogues. For 2tUra, experimental studies using transient absorption27−29 and time-resolved photoelectron spectroscopy30,23 reported subpicosecond ISC (τISC between 200 and 800 fs depending on excitation energy and environment31,32) from the initial bright state via a dark intermediate state. The lifetime of the formed triplet state is 50−200 ps in the gas phase and 70−300 ns in solution.31 These experiments were interpreted through static quantum chemistry calculations,33−36 nonadiabatic dynamics simulations,37,38 and the simulation of (time-resolved) photoelectron spectra.32,39 The computations showed that the relaxation of the initially excited 1πSπ* state (S2) first leads to the lower 1nSπ* state (S1), from where ISC to the 3πSπ* state occurs (T 1 ). The absence of accessible S 1 /S 0 conical intersections, the smallness of the singlet−triplet energy gaps, and the non-negligible spin−orbit couplings (SOCs) were advocated to explain the high ISC rate and quantum yield.22 The fact that the triplet lifetime is 1000 times longer in solution than
1. INTRODUCTION The genetic information in DNA and RNA is encoded by the sequence of the canonical nucleobases adenine, guanine, cytosine, thymine (only in DNA), and uracil (only in RNA). However, besides these canonical nucleobases, nature employs a large set of nucleobase analogues for various functions.1,2 Modified nucleobases can also be exploited as fluorescence markers,3,4 or as antiviral,5 anticancer,6,7 immuno-suppression,8 or other9 drugs. Their broad applicability of nucleobase analogues stems from their structural similarity to the canonical nucleobases10 in combination with some properties (e.g., photoreactivity) that are clearly different in analogues compared to the canonical bases. An important example of nucleobase modification is given by uracil and its substituted analogues 2-thiouracil (2tUra) and 2selenouracil (2seUra). Both 2tUra and 2seUra can be found in the tRNA11−13 of various microorganisms.14 Here, 2seUra mostly present as 5-methylaminomethyl-2seUraseems to be involved in codon−anticodon interactions15,16 and appears to help protecting tRNA from oxidative stress.17 Furthermore, both 2tUra and 2seUra exhibit thyroid peroxidase inhibitor activity, with 6-n-propyl-2tUra being a clinically approved drug9 and 6-n-propyl-2seUra showing an even higher activity.18 Furthermore, 2seUra can also effectively inhibit Type I iodothyronine deiodinase.19,20 Besides their biological applications, the photophysical and photochemical properties of 2tUra and 2seUra have attracted considerable attention. Thio-substituted nucleobases are known to exhibit a photoreactivity that is very different from their parent compounds.21,22 After excitation by UV radiation, the canonical nucleobases rapidly deactivate back to the electronic © 2019 American Chemical Society
Received: March 1, 2019 Published: April 30, 2019 3730
DOI: 10.1021/acs.jctc.9b00208 J. Chem. Theory Comput. 2019, 15, 3730−3742
Article
Journal of Chemical Theory and Computation in vacuum was rationalized by the presence of two T1 minima that are differently affected by the environment.31,36 The abundance of results for 2tUra is contrasted with the scarcity of studies on the photophysics of 2seUra. To the best of our knowledge, neither time-resolved spectroscopy nor theoretical calculations have been carried out on 2seUra to date. One computational study on substituted thymidines40 has reported vertical excitation energies and Franck−Condon SOC matrix elements for the related 2-selenothymidine but no relaxation pathways. Quite interesting is a recent transient absorption study on 6-selenoguanine41 which reports that seleno-substitution retains ultrafast ISC to the triplet states and accelerates ISC back to the ground state. The present work is the first computational study on the excited-state relaxation dynamics of 2seUra. We present both static calculationsi.e., vertical excitation computations, optimization of critical points, and potential energy scans and nonadiabatic dynamics simulations. For the sake of comparison to 2tUra, these calculations were carried out analogously to refs 35 and 38. Two fundamental questions are addressed in this work. First, how do the larger relativistic effects and the lower ionization potential of Se (compared to S) affect the potential energy surface landscape in 2seUra and its photodynamics? This question is of general interest to understand the relationship between molecular structure and excited-state dynamics. Second, which photoreactivity can be expected for 2seUra? This issue is important to judge whether 2seUra could be either used as a medicinal druge.g., in place of 2tUra derivativeswithout light-induced side effects or as a photosensitizer for photodynamic therapy.29
Figure 1. Structure of 2seUra with atom numbering and active space orbitals. The larger CAS(14,10) is denoted by a blue frame, while the smaller CAS(12,9) is shown with a yellow frame. The indices of the orbitals refer to the atom or part of the molecule where the largest part of the orbital is localized.
2. COMPUTATIONAL DETAILS In our investigations, we focus exclusively on the seleno-oxo tautomeric form of 2seUra (see Figure 1 top right corner). This tautomer was reported to be the most stable form of 2seUra in gas phase and in aqueous solution,42 much in the same way as the thio-oxo form is the most stable tautomer of 2tUra.43 The following subsections provide computational details of the electronic structure levels of theory, the static calculations, and the dynamics simulations. 2.1. Electronic Structure Methods. Two electronic structure methods were employed: multistate complete active space perturbation theory of second order (MS-CASPT2)44 and the algebraic diagrammatic construction scheme to second order for the polarization propagator (ADC(2)).45,46 2.1.1. MS-CASPT2. We employed MS-CASPT244 as the reference level of theory, as it has been shown before35,38 to be adequate for thio-nucleobases, and we expect the same for 2seUra. The MS-CASPT2 calculations are based on molecular orbitals from complete active space self-consistent field (CASSCF)47 calculations with an active space of 12 electrons in 9 orbitals (denoted CAS(12,9)). The active space orbitals are shown in Figure 1 and include the 8 orbitals of the π system and the lone pair of the selenium atoms. For the vertical excitation calculation (see below) we also included the lone pair of oxygen (CAS(14,10)), but this orbital was removed from the active space in all optimizations and scans because it did not affect the energies notably but introduced significant convergence problems. The employed basis set is of ANO-RCC48 type and uses a valence double-ζ plus polarization contraction (ANO-RCCVDZP). Scalar-relativistic one-electron integrals were computed
with the Douglas-Kroll-Hess formalism.49 As Molcas cannot state-average over states of different multiplicities, we employed two independently optimized active spaces, one optimized for 4 singlet states and another one for 3 triplet states. The same number of states was also used for the multistate CASPT2 treatment.44 The MS-CASPT2 computations were performed with an IPEA shift50 of zero51 and an imaginary level shift of 0.3 au in order to avoid intruder states.52 Additionally, Cholesky decomposition53 of the two-electron integrals was used to speed up the calculations. The MS-CASPT2(14,10) vertical excitation calculation included 9 singlets and 8 triplets and employed the larger ANO-RCC-VTZP basis set in combination with an IPEA shift of 0.25 au51 and an imaginary level shift of zero. All SA-CASSCF and MS-CASPT2 calculations were performed with Molcas 8.0.54 2.1.2. ADC(2). The second electronic structure method employed was ADC(2).45,46 This level of theory is computationally much cheaper than MS-CASPT2especially due to the availability of analytical gradientsbut still provides reasonably accurate potential energy surfaces, as is shown below. The computations employed the def2-TZVP55 basis set, except for the dynamics simulations that employed def2-SVP55 instead. Note that only nonrelativistic one-electron integrals were employed. The ADC(2) calculations, as well as the corresponding MP2 ground state computations, were carried out with Turbomole 7.0.56 As reported elsewhere,38 the ADC(2) implementation in Turbomole uses a very efficient resolution-ofidentity (RI) implementation to speed up the calculations,57 which is therefore also used here. The appropriateness of these settings is justified in Section 3. 3731
DOI: 10.1021/acs.jctc.9b00208 J. Chem. Theory Comput. 2019, 15, 3730−3742
Article
Journal of Chemical Theory and Computation
Table 1. Vertical Excitation Energies of 2seUra in eV (Oscillator Strength in Parentheses) Computed at the MS(9,8)CASPT2(14,10)/ANO-RCC-TZVP, MS(4,3)-CASPT2(12,9)/ANO-RCC-DZVP, and ADC(2)/def2-TVZP Levels of Theoryb 2-selenouracil
2-thiouracil
charactera
MS-CASPT2(14,10) E (fosc)
MS-CASPT2(12,9) E ( fosc)
ADC(2) E (fosc)
MS-CASPT2(12,9)35 E (fosc)
ADC(2)38 E ( fosc)
nSeπ*2 πSeπ*2 1 πSeπ6* 1 nSeπ6* 1 nOπ4* 1 πOπ4* 1 nOπ6* 1 πOπ6* 1 nSeπ*Ryd 3 πSeπ*2 3 nSeπ*2 3 πOπ*6 3 π5π6*
3.46 (0.00) 3.98 (0.62) 4.45 (0.03) 4.51 (0.00) 5.16 (0.00) 5.32 (0.04) 6.03 (0.00) 6.12 (0.34)
3.22 (0.00) 3.75 (0.50) 4.12 (0.14)
3.29 (0.00) 3.97 (0.38) 4.66 (0.03) 4.49 (0.00) 4.95 (0.00) 5.34 (0.04) 5.40 (0.00)
3.70 (0.00)
3.75 (0.00) 4.92 (0.11) 4.42 (0.38) 5.28 (0.00)
3.18 3.46 3.99
2.94 3.27
1 1
5.91 (0.02) 3.08 3.20 3.91
3.76
4.30 (0.11) 4.89 (0.00)
4.58 (0.00)
3.24 3.76 3.83
5.91 (0.02) 3.42 3.64 4.40 3.98
For the states of 2tUra, nSe and πSe are replaced by nS and πS. bFor comparison, corresponding values for 2tUra at the MS(3,3)-CASPT2(12,9)/ccpVDZ35 and ADC(2)/aug-cc-pVTZ levels of theory.38 a
2.2. Static Calculations. Static calculations were employed to explore the potential energy surfaces of 2seUra. We optimized the ground state (S0) minimum as well as several excited-state minima (S1, 2× S2, T1) and minimum energy crossing points (MXPs) (2× S1/S2, S1/T2, T1/T2, T1/S0). The optimized S0 minimum was used to compute vertical excitation energies. The optimizations were carried out with the external optimizer feature of ORCA 3.0.3,58 which was fed with the appropriate gradients from either ADC(2) or MS-CASPT2. For the MXP optimizations, the gradients to follow were obtained with either a projection algorithm59 for singlet−triplet crossings or a penalty function method60 for conical intersections. Initial geometries were prepared from the corresponding geometries of 2tUra.30,35 The optimizations resulted in ten critical points for ADC(2) and ten critical points for MS-CASPT2. The corresponding Cartesian coordinates are given in the Supporting Information. Subsequently, the optimized critical points were connected in an appropriate order by linear interpolation in internal coordinates (LIIC) scans. These scans were used to check for the presence or absence of barriers. Note that these scans do not provide minimum-energy paths but only upper limits for the involved barriers.61 2.3. Nonadiabatic Dynamics Simulations. In order to investigate the temporal evolution of 2seUra after excitation, we performed nonadiabatic dynamics simulations with the SHARC (surface hopping including arbitrary couplings) method.62,63 This method is an extension of the popular Tully surface hopping method64 able to treat ISC and is the same method that was previously used for the simulations of 2tUra dynamics.37,38 The dynamics simulations employed the RI-ADC(2)/def2SVP level of theory. To obtain initial conditions for the simulations, we sampled 1000 geometries from the Wigner distribution of the harmonic ground state potential, using frequencies and normal modes from an MP2/def2-SVP frequency calculation. For each of these 1000 geometries, 8 singlet states were computed with ADC(2) to generate an ensemble-broadened absorption spectrum. This spectrum was obtained as the sum over Gaussians (full width at half-maximum of 0.2 eV, the smallest value that sufficiently smoothed the spectrum) centered at the excitation energies of all 8000 states, with height proportional to the oscillator strength. The initial
excited state was selected stochastically in the excitation window of 3.95−4.05 eV (306−314 nm) based on the oscillator strength,65 i.e., brighter states were more likely to be selected as initial state. However, in order to reduce computation time and to avoid convergence problems from higher states, the dynamics simulations only included 3 singlet states (S0, S1, S2) and 2 triplet states (T1, T2). Hence, all 107 selected initial conditions started in the S2. The 107 trajectories were propagated with SHARC for a total simulation time of 700 fs, using a 0.5 fs nuclear propagation step. The electronic wave function was propagated with a 0.02 fs time step (using interpolated quantities) with the local diabatization formalism.66 Note that the electronic wave function was a linear combination of a total of 9 states (3 singlets, 2 × 3 triplets), as in SHARC the Ms components of the triplet states are treated separately. An energy-based decoherence scheme67 was applied.
3. RESULTS 3.1. Vertical Excitations. Table 1 presents an overview of the vertical excited states, together with previously published results of 2tUra.35,38 For 2seUra, three calculations were performed, one with MS(9,8)-CASPT2(14,10)/ANO-RCCTZVP, one with MS(4,3)-CASPT2(12,9)/ANO-RCC-VDZP, and one with ADC(2)/def2-TZVP. The ADC(2) calculation has the advantage that no bias is introduced by a limited active space. The large MS-CASPT2(14,10) calculation is taken as the reference, whereas the small MS-CASPT2(12,9) and ADC(2) settings are the levels of theory at which the optimizations and scans were carried out. For all three methods, the lowest excited state at the Franck− Condon geometry has 1nSeπ*2 character (see Figure 1 for the orbitals). This state is hence an excitation localized on the selenourea part of the molecule. At an energy of about 3.3 eV, this state is dark and hence does not contribute to the absorption spectrum. In this regard, the state is very similar to the 1nSπ*2 state of 2tUra, although the 1nSeπ2* is about 0.5 eV lower in energy. The second excited state is the 1πSeπ*2 state, with a vertical excitation energy of about 4.0 eV (3.8 eV with MS-CASPT2(12,9)) and an oscillator strength of 0.4−0.6. Again, this is a local excitation on the selenourea group. Interestingly, this state 3732
DOI: 10.1021/acs.jctc.9b00208 J. Chem. Theory Comput. 2019, 15, 3730−3742
Article
Journal of Chemical Theory and Computation
Figure 2. Prototypical geometries of the critical points of 2seUra, optimized with ADC(2)/def2-TZVP and MS-CASPT2(12,9)/ANO-RCC-VDZP. The large labels directly below the geometries indicate the main geometric features. The smaller, framed labels indicate which optimized critical points are similar to the shown prototypical geometries, e.g., the MS-CASPT2 S1 minimum or the T1/T2 MXP have geometries similar to (b2). The labels “ADC” and “PT2” indicate that a critical point is similar to different prototypes depending on the level of theory. See Tables 2 and 3 for corresponding geometrical parameters and the Supporting Information for Cartesian coordinates.
The ground state geometry of 2seUra is nearly planar with both the MP2 and MS-CASPT2 methods, as shown in Figure 2a. The aromatic ring has virtually the same bond lengths and angles as in 2tUra.35 However, the C−Se bond is significantly longer (1.79 Å) than the C−S bond in 2tUra (1.67 Å). In contrast to the ground state geometry, none of the excitedstate critical points are planar. Hence, we quantify the extent of ring deformationcalled ring puckeringby the Cremer-Pople parameter68 Q (amplitude) and the classification scheme of Boeyens.69 The latter scheme allows identifying the type of ring conformation, e.g., envelope (E), boat (B), screw-boat (S), or twist-boat (T). The indices given in the Boeyens symbols, e.g., B3,6, indicate which atoms (atom numbering in Figure 1) are displaced away from the ring plane. The displacement of the selenium atom from the ring plane is described by the pyramidalization angle of the Se-CN2 group, defined as 90° minus the angle between the C−Se bond and the normal vector of the CN2 plane.35 Both ADC(2) and MS-CASPT2 predict a minimum of the
seems to be in contrast to 2tUra, where the thiourea-localized ππ* state (1πSπ2*) is darker than the charge-transfer ππ* state (1πSπ6*).35 The corresponding state of 2seUrathe 1πSeπ6* stateshows only a smaller oscillator strength of 0.03−0.14 and is located at energies of 4.1−4.7 eV, depending on the method. Both states are located at lower energies than their 2tUra counterparts, indicating that the absorption spectrum of 2seUra is significantly red-shifted. The fourth excited singlet state is the 1nSeπ*6 state, which has an energy of 4.5 eV and is dark. Hence, the four lowest singlet states are all excitations from selenium orbitals. Excitations from other occupied orbitalslocalized on the oxygen or the ring only appear above 5 eV. Hence, for the present investigations we will focus only on the excitations from the selenium orbitals. The two lowest triplet states have 3πSeπ2* and 3nSeπ2* characters, making them the triplet analogues of the two lowest singlet states. Their energies are about 3.1−3.4 eV. Higher triplet states only appear after a notable energy gap, again allowing us to focus on the lowest triplet states. When comparing the three levels of theory, a generally good agreement is found. The order of the excited states is mostly reproduced, although ADC(2) reorders the close-lying 1nSeπ*6 and 1πSeπ6* states. Due to the smaller number of states included and the neglect of the nO orbital, the MS-CASPT2(12,9) calculation cannot describe most of the higher states. The MSCASPT2(12,9) also produces a red-shift of about 0.2 eV for all excited states, but the energy gaps between the states are in satisfactory agreement with the large MS-CASPT2(14,10) calculation. The ADC(2) and MS-CASPT2 results are also in good agreement with previously reported vertical excitation energies at the M06/6-31+G* level of theory for 2selenothymidine40 (S1: 3.43 eV, S2: 3.97 eV, T1: 2.91 eV, T2: 3.28 eV). We thus conclude that the computationally efficient ADC(2) and MS-CASPT2(12,9) levels of theory provide sufficiently accurate results compared to the MS-CASPT2(14,10) level of theory. 3.2. Excited State Minima and Crossing Points. The optimized critical points of 2seUra are depicted in Figure 2. Important geometry parameters and relative energies are given in Table 2 for ADC(2) and in Table 3 for MS-CASPT2.
nSeπ*2 state, which is the S1 both at the Franck−Condon point and at the minimum geometry. The most striking feature of this geometry is the large pyramidalization angle of the selenourea group of about 38°. This pyramidalization is mainly due to the electron in the π*2 orbital, which effectively changes the C2 atom from sp2 to sp3 hybridization and thus induces nonplanarity. The geometry differs slightly between the two methodswith ADC(2), the selenium is in the molecular plane and the C2 atom is puckered strongly (Figure 2b1), whereas with MSCASPT2, the ring is planar and the selenium is moved out of the ring plane to keep the pyramidalization angle unchanged (Figure 2b2). This is reminiscent of the situation in 2tUra,37 where MSCASPT2 predicts the planar ring plus displaced sulfur geometry at the minimum, but the C2-puckered ring with in-plane sulfur is only 0.03 eV higher in energy. Hence, the two different-looking S1 minima from ADC(2) and MS-CASPT2 are simply a consequence of the very flat S1 potential. Another interesting finding is that most bond lengths in the S1 minimum are very similar to the S0 minimum, indicating that no bond length 1
3733
DOI: 10.1021/acs.jctc.9b00208 J. Chem. Theory Comput. 2019, 15, 3730−3742
Article
Journal of Chemical Theory and Computation
Table 2. Energies, Geometry Parameters, and SOCs for Excited-State Minima and Crossing Points of 2seUra at the ADC(2)/ def2-TZVP Level of Theoryb ADC(2) Figure 2 E (eV) r1,2 (Å) r2,3 (Å) r3,4 (Å) r4,5 (Å) r5,6 (Å) r6,1 (Å) r2,7 (Å) a1,2,7 (deg) a3,4,8 (deg) p7,2,1,3 (deg) p12,6,5,1 (deg) Q (Å) Boeyens SOC (cm−1)
S0 min cs
S1 min 1 nSeπ2*
Spyr 2 min 1 πSeπ2*
Sboat mina 2 1 πSeπ6*
T1 min 3 πSeπ2*
1 S1/Spyr nSeπ2*/ 2 1 πSeπ2*
1 S1/Sboat nSeπ2*/ 2 1 πSeπ6*
S1/T2 1nSeπ2*/ 3 πSeπ2*
T2/T1 3nSeπ2*/ 3 πSeπ2*
T1/S0 3πSeπ2*/ cs
(a) 0.00 1.37 1.36 1.41 1.45 1.35 1.37 1.79 122.4 120.2 0.0
(b1) 2.86 1.40 1.41 1.40 1.45 1.36 1.36 1.91 116.4 120.5 37.6
(b2) 3.15 1.39 1.39 1.40 1.45 1.36 1.36 1.97 110.7 120.9 48.6
(c) 3.60 1.32 1.35 1.48 1.41 1.41 1.38 1.91 120.8 115.3 −5.3
(b2) 2.64 1.40 1.40 1.39 1.46 1.36 1.36 1.91 108.0 121.2 41.4
(b1) 3.19 1.38 1.38 1.41 1.45 1.36 1.36 1.99 107.2 121.2 56.2
(d) 3.74 1.30 1.30 1.52 1.39 1.43 1.46 1.92 117.1 114.2 −1.2
(b2) 2.90 1.40 1.39 1.40 1.46 1.35 1.37 1.90 111.7 120.8 31.4
(b2) 2.70 1.40 1.40 1.40 1.45 1.36 1.36 1.92 115.1 120.5 38.6
(e) 3.23 1.43 1.42 1.40 1.45 1.35 1.36 1.99 102.0 120.1 66.6
0.0
0.8
0.2
5.1
−0.5
0.8
32.1
−0.7
−0.3
−0.6
0.00 planar
0.26 E2
0.18 E2
0.08 B3,6
0.01 planar
0.20 E2
0.09 5 S6
0.17 3 S2 650
0.04 6 T2
0.08 E2 ≫0
The “Sboat min” is not a true minimum at the ADC(2) level of theory but is listed here for completeness. bQ is the puckering amplitude,68 and 2 “Boeyens” refers to the configuration of the ring69 (E: envelope, B: boat, S: screw-boat, T: twist-boat). Sub-/superscripts of r (bond length), a (bond angle), p (pyramidalization angle), and “Boeyens” refer to the atom numbering given in Figure 1. a
Table 3. Energies, Geometry Parameters, and SOCs for Excited-State Minima and Crossing Points of 2seUra at the MSCASPT2(12,9)/ANO-RCC-VDZP Level of Theorya MSCASPT2
S0 min cs
S1 min 1 nSeπ2*
Spyr 2 min 1 πSeπ2*
Sboat min 2 1 πSeπ6*
T1 min 3 πSeπ2*
1 S1/Spyr nSeπ2*/ 2 1 πSeπ2*
1 S1/Sboat nSeπ2*/ 2 1 πSeπ6*
S1/T2 1nSeπ2*/ 3 πSeπ2*
T2/T1 3nSeπ2*/ 3 πSeπ2*
T1/S0 3πSeπ2*/ cs
Figure 2 E (eV) r1,2 (Å) r2,3 (Å) r3,4 (Å) r4,5 (Å) r5,6 (Å) r6,1 (Å) r2,7 (Å) a1,2,7 (deg) a3,4,8 (deg) p7,2,1,3 (deg) p12,6,5,1 (deg) Q (Å) Boeyens SOC (cm−1)
(a) 0.00 1.38 1.37 1.41 1.45 1.36 1.38 1.82 122.4 119.8 −0.1
(b2) 2.88 1.40 1.40 1.40 1.46 1.36 1.37 1.99 115.4 120.4 38.0
(b2) 3.18 1.37 1.37 1.41 1.45 1.36 1.38 2.12 110.0 120.6 48.5
(c) 3.32 1.32 1.34 1.51 1.40 1.42 1.42 1.93 120.6 114.6 −1.8
(b2) 2.55 1.39 1.40 1.40 1.46 1.36 1.36 2.00 110.3 120.9 45.1
(b2) 3.23 1.36 1.36 1.42 1.45 1.36 1.38 2.20 106.5 120.6 55.7
(d) 3.32 1.33 1.32 1.52 1.40 1.42 1.45 1.92 119.1 114.6 0.0
(b2) 2.93 1.40 1.39 1.40 1.46 1.36 1.38 1.96 105.7 120.9 39.8
(b2) 2.67 1.39 1.38 1.40 1.46 1.36 1.38 2.05 114.6 120.7 40.9
(e) 2.78 1.43 1.40 1.39 1.46 1.35 1.36 1.99 94.7 120.7 58.1
0.0
−0.9
−0.2
−0.2
19.2
0.2
24.8
−1.1
−1.4
−1.1
0.03 planar
0.11 3 S2
0.05 E2
0.11 B3,6
0.05 1 S2
0.07 E2
0.06 E6
0.14 E2 650
0.16 E2
0.10 3 S2 660
Q is the puckering amplitude,68 and “Boeyens” refers to the configuration of the ring69 (E: envelope, B: boat, S: screw-boat, T: twist-boat). Sub-/ superscripts of r (bond length), a (bond angle), p (pyramidalization angle), and “Boeyens” refer to the atom numbering given in Figure 1. a
similar to the S 1 minimum, as it also shows strong pyramidalization of the selenourea group; thus, we denote this minimum as Spyr 2 minimum. The other geometry exhibits a slight boat conformation (B3,6 with puckering amplitude Q = 0.08 Å) and bond length alteration in the ring (Figure 2c). This minimum. geometry is termed Sboat 2 The last obtained excited-state minimum is located on the T1 potential energy surface, corresponding to 3πSeπ*2 character. It is very similar to the S2pyr minimum with regard to the pyramidalization angle, although the ring is more planar in the T1 minimum. Interestingly, no second T1 minimum was
alteration takes place. The only exception is the C−Se bond, which becomes 0.1−0.2 Å longer in the excited state. The MS-CASPT2 optimizations found two minima on the S2 adiabatic surface: one corresponding to 1πSeπ2* and the other one corresponding to 1πSeπ6*. At the ADC(2) level, only the 1πSeπ2* minimum is a true minimum, whereas the optimization of the 1 πSeπ6* energy follows a very small gradient until eventually converging to the 1πSeπ2* minimum (for a side-by-side discussion of both methods, Table 2 includes a geometry from the smallgradient region). Of these two points, the 1πSeπ2* minimum is 3734
DOI: 10.1021/acs.jctc.9b00208 J. Chem. Theory Comput. 2019, 15, 3730−3742
Article
Journal of Chemical Theory and Computation
Figure 3. LIIC scans along optimized minima and crossing points of 2seUra, obtained with (a) ADC(2)/def2-TZVP//ADC(2)/def2-TZVP or (b) MS-CASPT2(12,9)/ANO-RCC-VDZP//MS-CASPT2(12,9)/ANO-RCC-VDZP. The relaxation pathway is marked with open black circles, optimized minima with full black circles, and involved crossing points with crosses. The black arrows show the vertical excitation from the Franck− Condon geometry to the second excited singlet state. The character of the optimized states is given below each minimum. Relevant computed SOCs (in cm−1) are given at each singlet/triplet crossing point.
metries) shown in panel (a) and the MS-CASPT2//MSCASPT2 scans in panel (b). We stress that linear interpolation scans do not find the minimum-energy path between two points and hence might overestimate barriers. However, since in Figure 3 all path segments are monotonic, all barriers are given by the energies of the critical points themselves, with no additional barriers in between. In the figure, the critical points are organized in three photophysically reasonable pathways, called Paths I, II, and III. Paths I and II both start from the Franck−Condon point (i.e., the S0 minimum) in the bright S2 (1πSeπ*) state. Path I then pyr leads first to the Spyr 2 minimum, continues to the S1/S2 MXP, and relaxes to the S1 minimum. On the contrary, Path II leads minimum and then goes through the S1/Sboat first to the Sboat 2 2 MXP to also reach the S1 minimum. Hence, both paths describe the relaxation from the Franck−Condon region until the S1 minimum. The main difference is that in Path I, the selenourea group already pyramidalizes while still on the S2 potential energy surface, whereas in Path II pyramidalization only sets in after switching to the S1 state. This distinction can be seen in the lower halves of panels (a) and (b) of Figure 3, which plot the evolution of the pyramidalization angle along the relaxation pathways. This different evolution of the pyramidalization angle naturally affects the energetics of the two pathways. There are some relevant differences between the ADC(2) and MSCASPT2 results. In Path I, the ADC(2) energies of the involved critical points are 3.97, 3.15, 3.19, and 2.86 eV, giving an initial energy release of 0.82 eV and a barrier of 0.04 eV to reach the conical intersection. For MS-CASPT2, the energy release is 0.57 eV, and the barrier is 0.05 eV. For Path II, ADC(2) gives an energy release of 0.37 eV and a barrier of 0.14 eV; MS-CASPT2
identified for 2seUra in the present work, neither for ADC(2) nor for MS-CASPT2. Such a second T1 minimum was found previously for 2tUra,35,38 with a 3πSπ*6 wave function character and a B3,6 boatlike geometry (analogous to Figure 2c). For 2seUra, all attempts of optimizing such a minimum converged easily to the pyramidalized T1 minimum. We were also able to optimize five relevant crossing points of 2seUra, including two S1/S2 MXPs. The first S1/S2 MXP corresponds to the minimum of the crossing seam between 1 nSeπ2* and 1πSeπ2* and is termed S1/Spyr 2 MXP. Its geometry is very similar to the S1 minimum. The second MXP, called S1/Sboat 2 MXP, connects the 1nSeπ*2 and 1πSeπ*6 states. Its geometry is similar to the Sboat 2 , although with additional pyramidalization of the H−C6CN group, as depicted in Figure 2d. ISC from singlets to triplets might be facilitated by the presence of a crossing between 1nSeπ*2 and 3πSeπ*2 at a geometry very similar to the S1 minimum. At this geometry, the 1nSeπ*2 is S1 and the 3πSeπ2* is T2, hence we term this the S1/T2 MXP. Subsequent decay to the T1 surface is enabled by a T1/T2 crossing between the 3πSeπ2* and 3nSeπ2* states. Note, however, that this crossing is passed diabatically, since both the T1 minimum and the S1/T2 MXP involve the 3πSeπ*2 , whereas the 3 nSeπ2* state is involved in neither of the structures. The final optimized point is a T1/S0 MXP between the closedshell configuration and the 3πSeπ*2 state. The geometry is shown in Figure 2e, displaying the largest pyramidalization angles of all optimized structures. 3.3. Potential Energy Surface Scans. In order to connect the different optimized geometries and to verify the absence of potential energy barriers, we carried out linear interpolation scans. These scans are presented in Figure 3, with the ADC(2) potential energies (based on the ADC(2) optimized geo3735
DOI: 10.1021/acs.jctc.9b00208 J. Chem. Theory Comput. 2019, 15, 3730−3742
Article
Journal of Chemical Theory and Computation
deactivation. It is also neither possible to estimate how fast IC and ISC actually occur nor whether the optimizations found all relevant critical points. This information can be extracted from complementary nonadiabatic dynamics simulations, which are presented below. 3.4. Simulated Absorption Spectrum. As described above in Section 3.1, there are several excited states with significant oscillator strengths. In order to generate the initial conditions for the excited-state dynamics simulations, it is important to identify these bright excited states for each of the 1000 initial geometries. To this end, we simulated the UV absorption spectrum of 2seUra at the ADC(2)/def2-SVP level of theory that we also use for the dynamics simulations (Figure 4).
delivers values of 0.43 eV and a barrier of
