With this work, we present a protocol for the parameterization of a Linear Vibronic Coupling (LVC) Hamiltonian for quantum dynamics using highly accurate multiconfigurational electronic structure methods such as RASPT2/RASSCF, combined with a maximum-overlap diabatization technique. Our approach is fully portable and can be applied to many medium-size rigid molecules whose excited state dynamics requires a quantum description. We present our model and discuss the details of the electronic structure calculations needed for the parameterization, analyzing critical situations that could arise in the case of strongly interacting excited states. The protocol was applied to the simulation of the excited state dynamics of the pyrene molecule, starting from either the first or the second bright state (S2 or S5). The LVC model was benchmarked against state-of-the-art quantum mechanical calculations with optimizations and energy scans and turned out to be very accurate. The dynamics simulations, performed including all active normal coordinates with the multilayer multiconfigurational time-dependent Hartree method, show good agreement with the available experimental data, endorsing prediction of the excited state mechanism, especially for S5, whose ultrafast deactivation mechanism was not yet clearly understood.

Quantum dynamical (QD) simulations in polyatomic molecules are often run with reduced dimensionality models, generated predetermining the most important coordinates on the grounds of chemical intuition. This approach is advantageous since it strongly reduces the computational effort necessary to generate high-dimensionality potential energy surfaces (PESs) and to run QD in many dimensions with traditional methods.1 On the other side, recent methodological advances have made possible also the propagation of wavepackets (WP) in many dimensions, opening the route to a non-phenomenological description of decoherence and energy redistribution. The methods of reference in this field are probably the multiconfigurational time-dependent Hartree (MCTDH)2–4 and its multilayer (ML) extension (ML-MCTDH).5–8 They are extremely effective, even for nonadiabatic problems, especially if the coupled PESs have some simple functional form, like a low-order Taylor expansion in normal coordinates. Hamiltonians that use these simplified PESs are often referred to as model vibronic coupling Hamiltonians.9,10 They use a diabatic representation and quadratic expansions for the diagonal and off-diagonal PESs. If no other approximation is invoked, the above definition describes what is known as the quadratic vibronic coupling (QVC) Hamiltonian. However, it is usually further assumed that all diagonal PESs share the same normal modes and frequencies (usually taken all equal to the ones of the initial state before photo-excitation) and that off-diagonal terms are linear functions of the coordinates. These assumptions lead to the so-called linear vibronic coupling (LVC) model. LVC is the simplest Hamiltonian that can describe Conical Intersections (CoIs) and their multidimensional extensions (intersection seams), and in fact, it can be seen as a generalization to many states and modes of the two-state two-mode model adopted long-time ago to investigate the CoI problem.11 Model vibronic Hamiltonians have been quite successful to introduce the effect of interstate couplings in electronic spectra and to clarify the main features of a nonadiabatic dynamics around a CoI.9,10 Despite the “model” attribute, they can be adopted also for accurate descriptions of realistic problems, especially if the investigated molecules are rigid and/or the timescale of interest is very short (∼100 fs). As a matter of fact, in the last decade, they have been employed in the study of fast intersystem crossings in metal–organic complexes, 12–15 in ππ*/nπ* decays in nucleobases,16–20 and also to couple QD simulations with an explicit description of the environment.18,21 It is further worth noting that also popular models for excitonic problems essentially belong to the same family of Hamiltonians,22–26 but they (also) include off-diagonal constant terms so that CoIs cannot occur and the adiabatic PES are characterized by avoided crossings.

It became increasingly evident that even in the ultrafast regime, the QD can be drastically dependent on the parameters of the vibronic Hamiltonians, especially if the investigated system is characterized by several coupled quasi-degenerate states. In fact, the rate and yield of the predicted non-radiative processes can be totally different employing different Density Functional Theory (DFT) functionals19,20 or even different descriptions of the environment.18 These findings highlight the necessity to work out effective protocols to parameterize model Hamiltonians with electronic structure methods as accurate as possible.

We recently proposed a method based on a maximum-overlap diabatization to parametrize LVC Hamiltonians with time-dependent DFT (TD-DFT) calculations,19 which is very effective also for several excited states (10–20) and molecules with many degrees of freedom (100).27 From the point of view of electronic calculations, it only requires the ability to run single-point calculations and compute the overlap between electronic wavefunctions (WFs) at geometries displaced along the normal modes. Therefore, in principle, it is suitable for many electronic structure methods, and indeed, it is inspired by a procedure formerly proposed for configuration interaction (CI) WFs.28 Moreover, its computational cost is similar to that required to obtain the numerical gradients of all the involved states, and since the necessary calculations are embarrassingly parallel, even accurate methods can be adopted.

Multiconfigurational methods based on complete active space self-consistent field (CASSCF) and subsequent perturbative corrections (CASPT2) and their generalized extensions RASSCF and RASPT2 are, at the state of the art, among the most reliable electronic structure methods for computational photophysics and photochemistry. One of their major qualities is the capability to treat with similar accuracy states with different nature, including charge-transfer and double-excited states that challenge TD-DFT, provided the active space is properly selected. However, the dependence of the results on the active space composition, on the number of electronic states, and on the form of the zeroth order Hamiltonian makes LVC parameterization based on the CASSCF/CASPT2 protocol a rather intricate task. In particular, the formulation of the Fock operator in the construction of the zeroth order Hamiltonian has spawned several flavors of the perturbative correction, multi-state (MS),29 extended multi-state (XMS),30 and, more recently, extended dynamically weighted31 CASPT2, as well as the single-state single-reference and multi-state multi-reference variations of the MS-CASPT2.32 

In this contribution, we present, at the best of our knowledge, the first LVC Hamiltonian parameterized with (X)MS-RASPT2/RASSCF calculations for a medium-size molecule, such as pyrene. Pyrene is an interesting molecule that exhibits absorption bands of different bright states with a clear vibronic structure in the deep UV. Its photoinduced dynamics is characterized by the ultrafast internal conversion (IC) to the lowest dark excited state. While the IC process from the first bright excited state (320 nm) has been studied in detail both experimentally33–36 and theoretically,37,38 the IC process from the second excited state has been addressed only recently with transient absorption and bidimensional and photoelectron spectroscopy.36,39,40 Thanks to the unprecedented time-resolution (down to 6 fs), transient spectroscopy has allowed to resolve quantum beatings due to the motion of the vibrational WP in the excited state. Still, the picture of the IC mechanism from the second bright state is incomplete. Picchiotti et al.39 and Noble et al.40 recognized the involvement of intermediate dark states, but their role in the IC is not well understood yet.

We will study the decay dynamics of pyrene photoexcited to either its first or second bright states, adopting LVC Hamiltonians that fully account for the couplings of the lowest seven excited states and include all the active nuclear coordinates (49). We will evaluate the reliability of LVC PES by recomputing energies at relevant points of the dynamics, such as minima and energy-accessible CoIs. Moreover, we will investigate in depth the dependence of the QD results on different parameterizations of the Hamiltonian obtained with different active spaces and different implementations of the perturbative corrections. A parameterization of an LVC Hamiltonian is, actually, a much more stringent test of the stability of the computational protocol than the computation of the vertical excitations and/or of the numerical gradients, and we will analyze our results to enunciate few recommendations for future studies.

We consider a n dimensional diabatic basis, |d⟩ = (|d1⟩, |d2⟩, …, |dn⟩), and the following expression of the Hamiltonian:

H=iK+Viidia(q)|didi|+i,j>iVijdia(q)|didj|+|djdi|,
(1)

where q is the column vector of the ground state (GS) dimensionless normal coordinates. According to the Linear Vibronic Coupling (LVC) model, the kinetic (K) and potential (V) terms have the following form:

K=12pTΩp,
(2)
Viidia(q)=Ei0+λiiTq+12qTΩq,
(3)
Vijdia(q)=λijTq,
(4)

where Ω is the diagonal matrix of the GS normal-modes frequencies, p is the vector of the conjugated momenta, and T indicates the standard transpose operation for matrices. Therefore, the diagonal terms of the potential energy Viidia(q) are described in the harmonic approximation, and they share the same frequencies as the GS. The linear terms in the Hamiltonian represent the diabatic energy gradients λii and the inter-state diabatic couplings λij (ij).

The LVC Hamiltonian is parameterized by defining diabatic states |di⟩ to be coincident with the adiabatic reference states |ai⟩ at a reference geometry. We choose the GS minimum as reference. At displaced geometries, diabatic states are defined so to remain as similar as possible to the reference states |a(0)⟩. This idea was already proposed by Cimiraglia et al.28 for configuration-interaction WFs and then extended to TD-DFT by Neugebauer et al.41 and by some of us.19 More precisely, we follow the derivation presented in Ref. 19, and for each displaced geometry 0 + Δα (since now on Δα), we compute the adiabatic states |aα)⟩ and the matrix Sα) of their overlaps with |a(0)⟩,

Sij(Δα)=ai(0)|aj(Δα).
(5)

The transformation matrix D that defines the diabatic states at Δα,

|d=|a(Δα)D(Δα),
(6)

is then obtained as

D=ST(SST)12,
(7)

where for brevity the dependence on Δα is not explicitly reported. In Eq. (7), a Löwdin orthogonalization is used to account for the fact that the set of the computed adiabatic states at the displaced geometries is finite and therefore not complete.

At each displaced geometry the computed adiabatic energies form a diagonal matrix Vad(Δα)=diag(E1ad(Δα),E2ad(Δα),,Enad(Δα)) and the diabatic potential terms are simply

Vdia(Δα)=DT(Δα)Vad(Δα)D(Δα).
(8)

Therefore, the gradients λii and couplings parameters λij can be obtained from numerical differentiation with respect to each qα,

λij(α)=Vijdia(q)qαVijdia(Δα)Vijdia(Δα)2Δα.
(9)

In the following, the normal coordinates q and frequencies Ω were obtained at the second order perturbation theory level (MP2), whereas the energies Eiad(Δα) of the adiabatic states at each displaced geometry and their overlap S with the wave functions at the reference geometry were obtained at the RASSCF/RASPT2 level.

The vibronic wavefuction is defined in terms of the diabatic basis as |Ψ(q, t)⟩ = i|di⟩|Ψi(q, t)⟩, and the time evolution is computed by solving the time-dependent Schrödinger equation,

i|Ψi(q,t)t=H|Ψi(q,t).
(10)

In the following, we will investigate the time evolution of the population of the diabatic states. For state i at time t, it is simply Pi(t) = ⟨Ψi(q, t)|Ψi(q, t)⟩.

Pyrene is a highly symmetric molecule (D2h symmetry) with 26 atoms and 72 normal modes (see Tables S1–S3 in the supplementary material). For the parameterization of the LVC Hamiltonian, we have identified our diabatic states with the lowest seven excited adiabatic states at the S0 equilibrium geometry, belonging to four different irreducible representations: one state in Ag, two in B3u, one in B2u, and three in B1g. Then, we have displaced the atoms along each normal coordinate (obtained at MP2/ANO-L-VDZP level) both in the positive and in negative direction and calculated two main quantities: excitation energies and WF overlaps Siref|Sjdispl between all the eigenstates at the displaced and reference geometry (details on the WF overlap calculations at different geometries are given in the supplementary material). These data are then utilized to parameterize the LVC Hamiltonian according to Eqs. (7)(9). We note that while energy gradients are present only along symmetry conserving (Ag) modes, interstate couplings also exist along modes belonging to B1g, B2u, and B3u irreducible representations, which decrease the symmetry of the system, as indicated in Table I. 23 modes do not couple the electronic states of interest and are, therefore, excluded from the model. Our previous experience in the parameterization of the LVC Hamiltonian from TD-DFT indicates that a shift Δ = 0.1 in dimensionless coordinates guarantees accurate and robust results.19,42 Since diabatic states are built so to preserve at all geometries their electronic character, in the following, they will be named with the D2h symmetry labels of the adiabatic states they coincide with at the S0 minimum. Adiabatic states, on the contrary, will be denoted with the usual nomenclature Sx with x = 1, 2, …, 7 in order of increasing energy. It is worthy to remark that different diabatization techniques are actually possible.43 A strategy based on a one-shot computation of energy, gradients, and nonadiabatic coupling vectors with multireference Configuration Interaction Singles (CIS) and Configuration Interaction Singles and Doubles methods has been recently presented and implemented in Surface Hopping Including Arbitrary Couplings (SHARC) code.44 “Energy-based” methods, which rely only on energies and not on WFs, are also very attractive, and their simplicity makes them well suited to be applied also in combination with accurate and time-consuming electronic-structure methods such as CASSCF,45 Extended Multi-Configuration Quasi-Degenerate Perturbation Theory (XMCQDPT2),46 and Equation-of-Motion Coupled-Cluster Singles and Doubles (EOM-CCSD).47 Their implementation is very straightforward when each mode can only couple two states,46 while in the more general case, they require a fitting of the parameters, e.g., minimizing the root mean square deviation of the original ab initio and the model adiabatic PES at a representative number of points. The method we apply here is computationally demanding but is fully general. Moreover, being based on the overlaps of the WFs, it allows a direct and detailed control of the electronic character of the diabatic PESs.

TABLE I.

Coupling of the reference states along symmetry-breaking modes. Forbidden interactions in D2h symmetry are possible between states falling in the same irreducible representation of the lower point groups.

Irreducible representationPoint group atClassification of D2h states
of modesdisplaced geometriesinto new irreducible representations
B3u C2v A1: 1Ag, 1B3u, 2Ag, 2B3u B1: 1B2u, 1B1g, 2B1g, 3B1g 
B2u C2v A1: 1Ag, 1B2u, 2Ag B2: 1B3u, 1B1g, 2B3u, 2B1g, 3B1g 
B1g C2h Ag: 1Ag, 2Ag, 1B1g, 2B1g, 3B1g Bu: 1B3u, 1B2u, 2B3u 
Irreducible representationPoint group atClassification of D2h states
of modesdisplaced geometriesinto new irreducible representations
B3u C2v A1: 1Ag, 1B3u, 2Ag, 2B3u B1: 1B2u, 1B1g, 2B1g, 3B1g 
B2u C2v A1: 1Ag, 1B2u, 2Ag B2: 1B3u, 1B1g, 2B3u, 2B1g, 3B1g 
B1g C2h Ag: 1Ag, 2Ag, 1B1g, 2B1g, 3B1g Bu: 1B3u, 1B2u, 2B3u 

Electronic structure calculations with D2h and with reduced symmetry were performed at the RASPT2/RASSCF/ANO-L-VDZP level of theory. The calculations encompass the lowest eight roots of pyrene, which, due to the use of symmetry, fall in different irreducible representations. Three active spaces were used: a minimal one consisting of the frontier eight π and eight π* orbitals (full-π), with up to quadruple excitations [denoted as RAS(4, 80, 04, 8)] and two extended active spaces encompassing four and eight extra-valence virtual orbitals of π* character with a higher angular quantum number, denoted RAS(4, 80, 04, 12) and RAS(4, 80, 04, 16), respectively. The RASSCF scheme in which all molecular orbitals are put in RAS1 and RAS3 (leaving RAS2 empty) has been benchmarked previously, demonstrating the need of a high RAS1/RAS3 excitation level.48 The “empty RAS2” active space construction recipe has already shown to give accurate results for pyrene.39 We note that the extra-valence orbitals, despite bearing some resemblance to Rydberg orbitals, are not suitable for describing Rydberg states (not present among the states below 5 eV). Their only role is to capture more dynamic correlation at the RASSCF level, which has been shown to significantly improve the agreement with experimental data.49–51Figure 1 shows the active orbitals.

FIG. 1.

Active orbitals for pyrene in D2h symmetry, for each irreducible representation (top label; representations Ag, B1g, B2u, and B3u have no active orbitals). Bottom row (dark gray): π orbitals (RAS1), middle row (light gray) π* orbitals (RAS3), and top row (white): virtual orbitals with higher angular momentum (RAS3). The orbitals marked with * were excluded from the MS(8:12) and XMS(8:12) calculations.

FIG. 1.

Active orbitals for pyrene in D2h symmetry, for each irreducible representation (top label; representations Ag, B1g, B2u, and B3u have no active orbitals). Bottom row (dark gray): π orbitals (RAS1), middle row (light gray) π* orbitals (RAS3), and top row (white): virtual orbitals with higher angular momentum (RAS3). The orbitals marked with * were excluded from the MS(8:12) and XMS(8:12) calculations.

Close modal

In all calculations, on top of the RASSCF results, we have applied different types of perturbative corrections: either single state (SS), multi state (MS), or extended multi state (XMS) RASPT2, always using an imaginary shift of 0.2 a.u. and setting the IPEA shift to zero. For a more compact notation, each calculation will be labeled SS(n:m), MS(n:m), or XMS(n:m) depending on the type of perturbative correction, where n and m refer to the number of orbitals in RAS1 and RAS3, respectively. For calculations with D2h symmetry (at the reference geometry and along Ag modes), we rely on SS(8:16) energies which are virtually identical to MS results when the states are energetically separated and more accurate than XMS energies that rely on an average Fock operator. The only exception are the three close lying states belonging to the B1g irreducible representation for which MS(8:16) and XMS(8:16) energies were also evaluated. The SS(8:16) energies at the reference geometry were used as a uniform reference. For calculations with lower symmetry, we rely on (extended) multistate energies and WFs with reduced active space [i.e., (X)MS(8:12)] due to the interaction of near-degenerate states (forbidden at D2h symmetry) and the increase of computational effort. To allow for consistency, the change of energy along symmetry-reducing modes, evaluated at the (X)MS(8:12) level, was added to the reference SS(8:16) energies. The only exception are A1 states at geometries with C2v symmetry obtained by displacing along B3u modes, which were computed at the (X)MS(8:16) level as smaller active spaces were found to give nonphysically large interstate couplings. Overlaps were computed with the perturbatively modified WFs, obtained either at the (X)MS(8:12) or (X)MS(8:16) level. Further details on the calculations of the overlaps are given in Sec. III of the supplementary material. All the QM computations were performed with OpenMolcas,52,53 applying Cholesky decomposition.

ML-MCTDH wavepacket propagation2–8 was performed with the Quantics package.54,55 The method is also implemented in the original MCTDH code distributed upon request by Meyer and co-workers at Heidelberg University. The seven lowest energy excited states and the 49 (out of 72) normal coordinates with the appropriate symmetry to have non-vanishing couplings were included for all the LVC parametrized diabatic PESs. The dimension of the primitive basis set, the number of single particle functions, and the structure of the ML-MCTDH trees are shown in Sec. IV of the supplementary material for each type of calculation, together with some convergence tests (Fig. S9). We used a variable mean field (VMF) scheme with a fifth-order Runge–Kutta integrator of 10−7 accuracy threshold. The wavepackets were propagated for a total time of 2 ps.

The lowest seven excited states of pyrene belong to four irreducible representations (Table II). Among these states, we identify two optically bright states—1B2u with dominant configuration H(OMO) → L(UMO) and 2B3u with dominant configurations H − 1 → L + H → L + 1—as well as several dark states. Importantly, the lowest excited state is optically dark and, thus, responsible for the characteristic fluorescence of pyrene of hundreds of nanoseconds.59,60 We note the presence of a doubly excited state of Ag symmetry in the vicinity of the second bright state evidencing the need of multiconfigurational methods.

TABLE II.

Vertical excitation energies and transition dipole moment module (TDM) at the reference geometry for the first seven excited states of pyrene, obtained with the full-π active space (8:8) and with the extended active spaces (8:12) and (8:16). States are labeled according to the irreducible representations of the D2h point group. In the third column are reported the most relevant configuration state functions (CFSs) describing each state (see Fig. 1 for the representation of the involved orbitals). The last column reports the experimental adiabatic transition energies in the gas phase56,57 for bright states or of two-photon absorption experiments in apolar solvent58 for dark states. The (8:16) active space results are all reported relative to the SS(8:16) ground state value.

TDMEnergy (eV)Experimental
StateLabelCSFs(Debye)SS(8:8)SS(8:12)SS(8:16)MS(8:16)XMS(8:16)ΔE0-0 (eV)
S0 1Ag GS … 0.00 0.00 0.00 … … … 
S1 1B3u H → L + 1 0.00 3.23 3.22 3.23 … … 3.3657  
  H-1 → L        
S2 1B2u H → L 1.83 3.55 3.69 3.75 … … 3.8456  
S3 1B1g H → L + 2 0.00 4.11 4.13 4.16 4.00 4.10 4.1258  
S4 2Ag (H → L)2 0.00 4.30 4.35 4.32 … … 4.2958  
S5 2B3u H → L + 1 1.73 4.18 4.35 4.43 … … 4.6656  
  H-1 → L        
S6 2B1g H-2 → L 0.00 4.28 4.46 4.56 4.64 4.48 4.5458  
  H → L + 2        
S7 3B1g H-3 → L 0.00 4.73 4.77 4.82 4.89 4.85 4.9458  
  H → L + 3        
TDMEnergy (eV)Experimental
StateLabelCSFs(Debye)SS(8:8)SS(8:12)SS(8:16)MS(8:16)XMS(8:16)ΔE0-0 (eV)
S0 1Ag GS … 0.00 0.00 0.00 … … … 
S1 1B3u H → L + 1 0.00 3.23 3.22 3.23 … … 3.3657  
  H-1 → L        
S2 1B2u H → L 1.83 3.55 3.69 3.75 … … 3.8456  
S3 1B1g H → L + 2 0.00 4.11 4.13 4.16 4.00 4.10 4.1258  
S4 2Ag (H → L)2 0.00 4.30 4.35 4.32 … … 4.2958  
S5 2B3u H → L + 1 1.73 4.18 4.35 4.43 … … 4.6656  
  H-1 → L        
S6 2B1g H-2 → L 0.00 4.28 4.46 4.56 4.64 4.48 4.5458  
  H → L + 2        
S7 3B1g H-3 → L 0.00 4.73 4.77 4.82 4.89 4.85 4.9458  
  H → L + 3        

The vertical excitation energies at the reference geometry, obtained at different levels of theory, are reported in Table II. The full-π (8:8) active space shows both quantitative and qualitative differences with respect to the stronger correlated (8:12) and (8:16) active spaces. Indeed, while the energies of states such as 2Ag, 1B1g, and 1B3u are already converged with respect to the active space size, the remaining states (in particular, both bright states 1B2u and 2B3u) exhibit strong dependence on the active space size, being red-shifted by 0.2 eV–0.3 eV at the SS(8:8) level with respect to SS(8:16). As a consequence of the unbalanced description, the energy order of the states changes as a function of the active space (Table II) with profound consequences for the QD simulations. The trend in the (8:8)-(8:12)-(8:16) sequence evidences that energies are not fully converged even with the largest active space, but they show an asymptotic behavior. Accordingly, comparison with the experimental gas-phase data56–58 shows that the computed transition energies of the bright states are underestimated. The SS(8:16) set provides closest agreement, thus implicitly supporting the predicted state order.

Concerning the type of perturbative correction, the SS-variation of the RASPT2 method is the best approximation with D2h symmetry where states of the same irreducible representation are far apart in energy and do not mix. Only in the case of the B1g irreducible representation, (X)MS-RASPT2 energies were considered due to the proximity of the electronic states. Indeed, the three methods predict energies that deviate by up to 0.16 eV. XMS-RASPT2, whose use is advocated for near-degenerate and strongly interacting electronic states,61 is found to deviate only marginally from the SS-RASPT2 results. Eventually, considering the computational cost and the small error, SS(8:16) was used to calculate the energies along symmetry-conserving normal modes.

At the S0 equilibrium geometry, all the excited states show a gradient only along the totally symmetric Ag modes. With the numerical gradients at hand, within the displaced harmonic oscillator approximation, we can predict the structures of the minima of the adiabatic states and the reorganization energies λ (details in the supplementary material). Interestingly, we obtain small reorganization energies (up to ∼0.3 eV, Table III), which reflects the rigidity of the pyrene molecule and justifies the harmonic approximation underlying the LVC model. The predicted structures and reorganization energies are in a very good agreement with the results from explicit optimizations at the SS-RASPT2/RASSCF(4, 80, 04, 8)/ANO-L-VDZP level39 [i.e., SS(8:8), Table III].62 Taking into consideration the reorganization energies resolves the apparent disagreement between experiment and theory regarding the energetic order of 2B3u and 2B1g (Table II). Two-photon absorption experiments put the 2B1g (4.54 eV) below the second bright state 2B3u (4.66 eV) at the respective excited minimum. When the reorganization energies—predicted as ∼0.05 eV for 2B3u and 0.23 eV for 2B1g (Table III)—are considered, the state order is inverted in the Franck–Condon (FC) point.

TABLE III.

Comparison between SS-RASPT2/RASSCF(4, 80, 04, 8)/ANO-L-VDZP optimized minima (OPT) and LVC model minima for the adiabatic excited states of pyrene: reorganization energy λ for each structure and RMSD between the two Cartesian structures for each state. The reorganization energies were obtained as the difference in energy between the reference geometry and the corresponding minimum at the SS(8:16) level (OPT) or by projecting the SS(8:16) gradient onto the normal modes (LVC, see the supplementary material).

S1S2S3S4S5S6
OPT39 LVCOPT39 LVCOPTLVCOPTLVCOPT39 LVCOPTLVC
λ (eV) 0.08 0.09 0.10 0.10 0.19 0.16 0.26 0.18 0.05 0.06 0.22 0.24 
RMSD 0.005 0.005 0.004 0.012 0.005 0.010 
S1S2S3S4S5S6
OPT39 LVCOPT39 LVCOPTLVCOPTLVCOPT39 LVCOPTLVC
λ (eV) 0.08 0.09 0.10 0.10 0.19 0.16 0.26 0.18 0.05 0.06 0.22 0.24 
RMSD 0.005 0.005 0.004 0.012 0.005 0.010 

Vibronic coupling between the considered diabatic states is observed both along totally symmetric Ag modes and along the symmetry-decreasing modes belonging to the B3u, B2u, and B1g irreducible representations (Table I). As noted earlier, in D2h symmetry, electronic states of the same irreducible representation are energetically well separated, which results in a weak interaction (coupling). On the other hand, displacement along symmetry-lowering modes allows also for interactions that were forbidden in D2h symmetry: this is particularly evident in the case of the first bright state S2, which is the only B2u state in D2h symmetry and otherwise would never be depopulated. Symmetry-lowering results in variable grouping of the states in irreducible representations of lower point groups. This requires a different state averaging along each of the three symmetry-decreasing sets of normal modes, which affects both the RASSCF and RASPT2 results, in particular, in the case of XMS-RASPT2, which relies on an average Fock operator. Moreover, the presence of close lying states requires the use of (X)MS-RASPT2 corrections. Because of this, the level of theory of the WF overlap calculations must be accurately selected for each irreducible representation of each point group so as to balance between computational cost and accuracy of the description. To assess the reliability of the reduced symmetry calculations in reproducing the electronic structure with the same precision as the D2h calculations, the electronic structure at the reference geometry was computed with each of the lower symmetries. Table IV and Fig. 2 show the deviation of the adiabatic energies at the (X)MS(8:12) and (X)MS(8:16) levels from the reference D2h-SS(8:16) values when the symmetry is reduced. The agreement with the reference values is generally good, with XMS-RASPT2 being more accurate than MS-RASPT2, which tends to overestimate the energy splitting and WF mixing in the case of strongly interacting states. Comparing the two active spaces, it is evident how the energies are sensitive to the degree of electronic correlation, with the (8:16) results being more faithful to the reference energies than the (8:12) ones, both for MS- and XMS-RASPT2. Thus, it is obvious that the best choice would be to calculate all the WF overlaps (necessary for the LVC parameterization) with the larger active space, but this is computationally very demanding. To balance between computational cost and accuracy of the description, we have computed the wavefunction overlaps at the (X)MS(8:12) level, except for critical situations (i.e., strongly interacting states), where we have used (X)MS(8:16), and that will now be discussed.

TABLE IV.

Vertical excitations at the reference geometry: deviation from the reference D2h-SS(8:16) values (reported in the first row) at different levels of theory. Positive and negative deviations larger than 0.10 in the absolute value are highlighted in bold and italic, respectively. For each symmetry, states of the same irreducible representation fall into the same RASPT2/RASSCF calculation. C2v(1) and C2v(2) refer to the reduced symmetry along modes B3u and B2u, respectively.

Deviation from reference energy (eV)Absolute meanStandard
SymmetryLevel of theoryS1S2S3S4S5S6S7deviation (eV)deviation (eV)
D2h SS(8:16) 3.23 3.75 4.16 4.32 4.43 4.56 4.82 … … 
Irreducible rep. A1 B1 B1 A1 A1 B1 B1   
C2v (1) MS(8:12) 0.04 −0.05 0.19 0.06 −0.06 0.08 0.09 0.082 0.12 
 MS(8:16) −0.01 −0.02 0.17 0.02 −0.03 0.10 0.08 0.061 0.10 
 XMS(8:12) 0.05 0.11 −0.04 0.05 −0.05 −0.04 0.11 0.064 0.09 
 XMS(8:16) 0.00 −0.08 −0.03 0.03 0.00 −0.02 0.09 0.036 0.06 
Irreducible rep. B2 A1 B2 A1 B2 B2 B2   
C2v (2) MS(8:12) 0.02 0.02 0.23 0.04 −0.05 0.10 0.10 0.080 0.13 
 MS(8:16) −0.01 0.02 0.20 0.01 −0.02 0.10 0.08 0.063 0.11 
 XMS(8:12) 0.01 0.00 −0.08 0.12 0.06 −0.06 0.09 0.060 0.09 
 XMS(8:16) −0.02 0.01 −0.06 0.09 −0.03 −0.06 0.06 0.047 0.07 
Irreducible rep. Bu Bu Ag Ag Bu Ag Ag   
C2h MS(8:12) 0.06 −0.06 0.22 0.03 −0.02 0.11 0.10 0.086 0.14 
 MS(8:16) 0.02 0.02 0.19 −0.03 −0.01 0.11 0.07 0.064 0.11 
 XMS(8:12) 0.07 0.10 −0.06 0.02 0.00 −0.04 0.09 0.055 0.08 
 XMS(8:16) 0.03 −0.03 −0.05 0.06 0.00 −0.06 0.04 0.039 0.06 
Deviation from reference energy (eV)Absolute meanStandard
SymmetryLevel of theoryS1S2S3S4S5S6S7deviation (eV)deviation (eV)
D2h SS(8:16) 3.23 3.75 4.16 4.32 4.43 4.56 4.82 … … 
Irreducible rep. A1 B1 B1 A1 A1 B1 B1   
C2v (1) MS(8:12) 0.04 −0.05 0.19 0.06 −0.06 0.08 0.09 0.082 0.12 
 MS(8:16) −0.01 −0.02 0.17 0.02 −0.03 0.10 0.08 0.061 0.10 
 XMS(8:12) 0.05 0.11 −0.04 0.05 −0.05 −0.04 0.11 0.064 0.09 
 XMS(8:16) 0.00 −0.08 −0.03 0.03 0.00 −0.02 0.09 0.036 0.06 
Irreducible rep. B2 A1 B2 A1 B2 B2 B2   
C2v (2) MS(8:12) 0.02 0.02 0.23 0.04 −0.05 0.10 0.10 0.080 0.13 
 MS(8:16) −0.01 0.02 0.20 0.01 −0.02 0.10 0.08 0.063 0.11 
 XMS(8:12) 0.01 0.00 −0.08 0.12 0.06 −0.06 0.09 0.060 0.09 
 XMS(8:16) −0.02 0.01 −0.06 0.09 −0.03 −0.06 0.06 0.047 0.07 
Irreducible rep. Bu Bu Ag Ag Bu Ag Ag   
C2h MS(8:12) 0.06 −0.06 0.22 0.03 −0.02 0.11 0.10 0.086 0.14 
 MS(8:16) 0.02 0.02 0.19 −0.03 −0.01 0.11 0.07 0.064 0.11 
 XMS(8:12) 0.07 0.10 −0.06 0.02 0.00 −0.04 0.09 0.055 0.08 
 XMS(8:16) 0.03 −0.03 −0.05 0.06 0.00 −0.06 0.04 0.039 0.06 
FIG. 2.

Vertical excitation energies at the reference geometry calculated with the reduced symmetries of the B3u modes (top right), B2u modes (bottom left), and B1g modes (bottom right). In the top left panel are reported the reference D2h-SS(8:16) energies. Full circles = S0 and bright states; empty circles = dark states. Vertical dotted lines connect states of the same irreducible representation for each point group and level of theory. The horizontal full lines set the reference D2h-SS(8:16) energies.

FIG. 2.

Vertical excitation energies at the reference geometry calculated with the reduced symmetries of the B3u modes (top right), B2u modes (bottom left), and B1g modes (bottom right). In the top left panel are reported the reference D2h-SS(8:16) energies. Full circles = S0 and bright states; empty circles = dark states. Vertical dotted lines connect states of the same irreducible representation for each point group and level of theory. The horizontal full lines set the reference D2h-SS(8:16) energies.

Close modal

For each group of symmetry-reducing modes, we can identify a pair of close lying states, which require particular attention, to make sure that the new state averaging scheme retains the relative state order and energy gaps as at the reference D2h geometry: S4/S5 along B3u modes [ΔESS(8:16)D2h = 0.11 eV], S5/S6 along B2u modes [ΔESS(8:16)D2h = 0.13 eV], and S3/S4 along B1g modes [ΔESS(8:16)D2h = 0.16 eV]. Table V shows the average, maximum, and minimum WF overlap (absolute value) for each critical couple of states. For S6–S5 (along B2u modes) and S4–S3 (along B1g modes), the (8:12) energy splitting is always overestimated with respect to the reference one, and the WF overlaps are consequently small, with XMS-RASPT2 being more accurate than MS-RASPT2. Even though from the theoretical point of view the overestimation of the energy gap is conceptually as wrong as its underestimation, from the practical point of view, a larger energy gap (which results in a smaller diabatic coupling in the final Hamiltonian) is not as dramatic as a too small energy gap since artificially large diabatic couplings can make the QD calculations much more problematic. On the contrary, the case of S4 and S5 states along B3u modes (i.e., A1 representation, see Fig. 2) is more critical: (X)MS(8:12) reduces the energy gap until near-degeneracy of the two states, producing an unphysically high WF overlap (and diabatic coupling, see Fig. S3 in the supplementary material for the correlation between accuracy of the ΔE and wavefunction mixing). Table V shows that at the MS(8:12) level, they are perfectly degenerate, resulting in an average WF overlap of about 0.40. On the other hand, increasing the active space, the energy gap increases, getting closer to the reference D2h-SS(8:16) value, and the S5–S4 mixing is significantly reduced [0.012 at MS(8:16) and 0.006 at XMS(8:16) level, see Table V].

TABLE V.

Energy gap and WF overlaps along symmetry reducing modes (average absolute value, minimum and maximum absolute values) between states S5–S4 (top), S6–S5 (middle), and S4–S3 (bottom) calculated with different symmetry and level of theory.

Deviation fromSiref|Sjdispl
ModesSymmetryLevel of theoryΔE (eV)reference ΔE (eV)AverageMinMax
S5S4 B3u C2v(1) MS(8:12) 0.00 −0.11 0.395 0.137 0.613 
   MS(8:16) 0.09 −0.02 0.012 0.001 0.044 
   XMS(8:12) 0.01 −0.10 0.080 0.001 0.262 
   XMS(8:16) 0.09 −0.02 0.006 8 × 10−5 0.020 
S6S5 B2u C2v(2) MS(8:12) 0.27 0.14 0.025 3 × 10−4 0.070 
   XMS(8:12) 0.13 0.00 0.029 0.001 0.112 
S4S3 B1g C2h MS(8:12) 0.41 0.25 0.030 0.005 0.090 
   XMS(8:12) 0.25 0.09 0.010 0.001 0.033 
Deviation fromSiref|Sjdispl
ModesSymmetryLevel of theoryΔE (eV)reference ΔE (eV)AverageMinMax
S5S4 B3u C2v(1) MS(8:12) 0.00 −0.11 0.395 0.137 0.613 
   MS(8:16) 0.09 −0.02 0.012 0.001 0.044 
   XMS(8:12) 0.01 −0.10 0.080 0.001 0.262 
   XMS(8:16) 0.09 −0.02 0.006 8 × 10−5 0.020 
S6S5 B2u C2v(2) MS(8:12) 0.27 0.14 0.025 3 × 10−4 0.070 
   XMS(8:12) 0.13 0.00 0.029 0.001 0.112 
S4S3 B1g C2h MS(8:12) 0.41 0.25 0.030 0.005 0.090 
   XMS(8:12) 0.25 0.09 0.010 0.001 0.033 

In conclusion, the (X)MS(8:12) WF overlaps represent a fair compromise between computational time and accuracy, except for the states of A1 representation along B3u modes (C2v symmetry), for which the bigger active space is needed to avoid artificially high S5/S4 overlaps. For comparison, we have produced three sets of data for the LVC parameterization: one in which all the overlaps were computed at the XMS(8:12) level and two sets in which the B3u-A1 states were computed with the bigger active space [i.e., MS(8:16) or XMS(8:16)].

The three different parameterizations of the LVC Hamiltonian will be named from now on LVCMS(16), LVCXMS(12), and LVCXMS(16) depending on the highest level of theory employed for the computation of the WF overlaps [MS(8:16), XMS(8:12) or XMS(8:16), respectively]. Figure 3 compares scans of the LVCMS(16) diabatic PESs along Ag collective coordinates leading from the S0 minimum to the minima of the different LVC diabatic PESs (solid lines) with the energies of the corresponding adiabatic states recomputed at the D2h-SS(8:16) level (scattered points). The comparison shows that LVC PESs are remarkably accurate, especially for the lower energy states. Some inaccuracies arise for 3B1g and 2B3u along the coordinate connecting the S0 and the 1B1g minima (Fig. 3, middle left panel). This is connected with the degeneracy, at distorted geometries, with a higher lying “intruder” state at the RASSCF level that is influencing the CASPT2 correction. We emphasize that upon (X)MS-CASPT2 correction, the “intruder” states blue-shift above 5 eV, which evidences that their involvement at the RASSCF level is merely an artifact of the unbalanced description of the electronic states when dynamic correlation is not considered.

FIG. 3.

Scans of the LVCMS(16) diabatic potential energy surfaces (dashed lines) along collective Ag coordinates connecting the 1Ag equilibrium geometry with the minima of the LVC diabatic states and corresponding adiabatic energies computed at the SS(8:16) level (hollow circles). Note that although the SS(8:16) states are adiabatic, they are distinguished by symmetry, which explains the observed crossings and justifies that for each symmetry, LVC adiabatic energies are very similar to LVC diabatic ones.

FIG. 3.

Scans of the LVCMS(16) diabatic potential energy surfaces (dashed lines) along collective Ag coordinates connecting the 1Ag equilibrium geometry with the minima of the LVC diabatic states and corresponding adiabatic energies computed at the SS(8:16) level (hollow circles). Note that although the SS(8:16) states are adiabatic, they are distinguished by symmetry, which explains the observed crossings and justifies that for each symmetry, LVC adiabatic energies are very similar to LVC diabatic ones.

Close modal

To have a closer look at the performance of the LVC model in the minima, we consider the LVCMS(16) parameterization and recomputed the SS(8:16) energies at all the diabatic minima located with the LVC model. Data in Table S8 of the supplementary material show that LVC and RASPT2 energies are extremely similar. The largest differences for a state in its own minimum are seen for 2Ag and 2B1g and are 0.04 eV. At each minimum, also the energies of the other states are quite similar with the partial exceptions of states 2B3u and 3B1g, which, far from their own minimum, can show an interaction with higher lying states at the RASSCF level not included in the model, as mentioned previously.

With the LVC model it is also possible to analytically determine the lowest energy crossing of pairs of diabatic states in D2h symmetry. Note that since in D2h off-diagonal couplings among states of the same symmetry are possible, diabatic and adiabatic LVC states do not coincide, and therefore, these crossings do not correspond, rigorously speaking, to CoIs between adiabatic states. However, we already showed that mixings between states of the same symmetry are minimal when the D2h point group is applied. Table VI reports the LVC and SS(8:16) energies of all states at crossings with energies lower than 4.5 eV (i.e., accessible from 2B3u, whose vertical excitation energy is 4.43 eV). For crossing up to 4.5 eV, the agreement is remarkably good. RASPT2 confirms that these structures correspond to points of quasi-degeneracy, and in most of the cases, also the LVC absolute energy is correct up to few hundredths of eV. In particular, LVC correctly predicts that the 1B1g/2B3u crossing actually corresponds to a quasi-triple CoI involving also the 2Ag state and reproduces the absolute energies up to 0.02 eV. A further quasi-triple CoI involving the 1B3u, 1B2u, and 1B1g states (proposed previously based on orbital analysis and CoI search39) is also confirmed. In this case, however, LVC overestimates the energy by ∼0.10 eV–0.15 eV. Considering diabatic crossings at higher energy (check Table S9 in the supplementary material), LVC predictions are still rather reliable, but, as expected, differences with respect to RASPT2 energies increase. Interestingly, LVC correctly predicts that at 1B3u/2Ag crossing, four states are found in <0.17 eV (i.e., also 1B2u and 1B1g), suggesting that a quasi-fourfold CoI might exist in the proximity of that structure.

TABLE VI.

Diabatic [LVCMS(16)] and adiabatic [RASPT2, SS(8:8)] energies (eV) of pyrene at a number of crossing points between LVC diabatic states. Bold characters highlight states that are quasi-degenerate (data for higher energy crossings are in Table S9 of the supplementary material).

States
S1S2S3S4S5S6S7
CoIMethods1B3u1B2u1B1g2Ag2B3u2B1g3B1g
1B3u/1B2u LVC 4.20 4.20 4.42 4.53 5.37 4.76 5.64 
 RASPT2 4.16 4.16 4.37 4.37 5.72 5.07 5.07 
1B3u/1B1g LVC 4.43 4.46 4.43 4.81 5.57 4.84 5.73 
 RASPT2 4.27 4.29 4.33 4.60 5.71 4.81 5.41 
1B2u/1B1g LVC 4.04 4.25 4.25 4.68 5.18 4.68 5.37 
 RASPT2 3.88 4.12 4.18 4.56 5.57 4.63 5.14 
1B1g/2B3u LVC 3.20 3.89 4.45 4.50 4.45 4.80 4.91 
 RASPT2 3.21 3.89 4.47 4.49 4.46 4.88 4.96 
2Ag/2B3u LVC 3.17 3.80 4.27 4.40 4.40 4.64 4.83 
 RASPT2 3.17 3.80 4.27 4.39 4.40 4.65 4.82 
2B3u/2B1g LVC 3.19 3.67 4.06 4.20 4.40 4.40 4.79 
 RASPT2 3.19 3.68 4.06 4.20 4.41 4.42 4.75 
States
S1S2S3S4S5S6S7
CoIMethods1B3u1B2u1B1g2Ag2B3u2B1g3B1g
1B3u/1B2u LVC 4.20 4.20 4.42 4.53 5.37 4.76 5.64 
 RASPT2 4.16 4.16 4.37 4.37 5.72 5.07 5.07 
1B3u/1B1g LVC 4.43 4.46 4.43 4.81 5.57 4.84 5.73 
 RASPT2 4.27 4.29 4.33 4.60 5.71 4.81 5.41 
1B2u/1B1g LVC 4.04 4.25 4.25 4.68 5.18 4.68 5.37 
 RASPT2 3.88 4.12 4.18 4.56 5.57 4.63 5.14 
1B1g/2B3u LVC 3.20 3.89 4.45 4.50 4.45 4.80 4.91 
 RASPT2 3.21 3.89 4.47 4.49 4.46 4.88 4.96 
2Ag/2B3u LVC 3.17 3.80 4.27 4.40 4.40 4.64 4.83 
 RASPT2 3.17 3.80 4.27 4.39 4.40 4.65 4.82 
2B3u/2B1g LVC 3.19 3.67 4.06 4.20 4.40 4.40 4.79 
 RASPT2 3.19 3.68 4.06 4.20 4.41 4.42 4.75 

Figure 4 shows the time evolution of the electronic populations up to 2 ps after the initial photo-excitation to either the first (1B2u) or the second (2B3u) bright states according to the LVCMS(16) and LVCXMS(16) parameterizations (results with LVCXMS(12) are given in Figure S14 of the supplementary material). The insets report a close-up of the same data in the first 100 fs. LVCMS(16) and LVCXMS(16) Hamiltonians deliver similar predictions: 1B2u decays essentially on the lowest state 1B3u, while after an initial excitation to 2B3u, we observe a fast (<20 fs) rise of a transient population of some intermediate states, followed by a only slightly slower population of the first bright state 1B2u, which reaches its maximum population (∼0.5) in 100 fs and then slowly decays toward 1B3u. The intermediate population of 1B2u is consistent with the two-step interpretation of Borrego-Varillas et al. who reported transient signatures of 1B2u when pumping the second bright state.36 Moreover, the delayed decay to the lowest excited state (on a 0.5 ps time scale) observed after excitation to 2B3u agrees with experimental time constants reported in the literature.36,39,40

FIG. 4.

Dynamics of the populations of the diabatic electronic states obtained by initially exciting the wavepacket on 1B2u (left) or 2B3u (right) states for the LVCMS(16) [panels (a) and (b)] and LVCXMS(16) [panels (c) and (d)] parameterizations. The insets highlight the dynamics in the first 100 fs.

FIG. 4.

Dynamics of the populations of the diabatic electronic states obtained by initially exciting the wavepacket on 1B2u (left) or 2B3u (right) states for the LVCMS(16) [panels (a) and (b)] and LVCXMS(16) [panels (c) and (d)] parameterizations. The insets highlight the dynamics in the first 100 fs.

Close modal

A closer analysis highlights some differences. For an excitation to 1B2u, the decay to 1B3u is faster according to LVCMS(16) than according to LVCXMS(16). Thereby, the LVCMS(16) dynamics agrees better with experiments, uniformly assigning a sub-100 fs time constant to the S2 → S1 IC. Analysis of the couplings (Table VII) suggests that this finding can partially arise from the larger coupling predicted by LVCMS(16) (norm: 0.043 eV) than by LVCXMS(16) (norm: 0.030 eV) and mainly due to the contribution of mode 60: 0.025 eV in LVCMS(16) and 0.010 eV in LVCXMS(16). However, further motivations will be highlighted below.

TABLE VII.

Norm of the diabatic coupling vectors for MS(8:16) and XMS(8:16) parameterizations. Bold numbers highlight differences between the two parameterizations that have a remarkable impact on the population dynamics.

MS(8:16)XMS(8:16)
State1B3u1B2u1B1g2Ag2B3u2B1g3B1g1B3u1B2u1B1g2Ag2B3u2B1g3B1g
1B3u 0.159       0.159       
1B2u 0.043 0.184      0.030 0.184      
1B1g 0.199 0.196 0.257     0.111 0.096 0.257     
2Ag 0.108 0.116 0.049 0.257    0.105 0.043 0.058 0.256    
2B3u 0.108 0.124 0.027 0.054 0.126   0.072 0.176 0.096 0.028 0.126   
2B1g 0.087 0.096 0.037 0.235 0.152 0.266  0.056 0.126 0.059 0.109 0.146 0.267  
3B1g 0.175 0.238 0.042 0.089 0.077 0.073 0.143 0.105 0.028 0.089 0.093 0.046 0.037 0.142 
MS(8:16)XMS(8:16)
State1B3u1B2u1B1g2Ag2B3u2B1g3B1g1B3u1B2u1B1g2Ag2B3u2B1g3B1g
1B3u 0.159       0.159       
1B2u 0.043 0.184      0.030 0.184      
1B1g 0.199 0.196 0.257     0.111 0.096 0.257     
2Ag 0.108 0.116 0.049 0.257    0.105 0.043 0.058 0.256    
2B3u 0.108 0.124 0.027 0.054 0.126   0.072 0.176 0.096 0.028 0.126   
2B1g 0.087 0.096 0.037 0.235 0.152 0.266  0.056 0.126 0.059 0.109 0.146 0.267  
3B1g 0.175 0.238 0.042 0.089 0.077 0.073 0.143 0.105 0.028 0.089 0.093 0.046 0.037 0.142 

For an excitation to 2B3u, the initial decay (∼10 fs) is toward 2B1g and 2Ag according to LVCMS(16) and toward 2B1g, 1B1g, and directly 1B2u according to LVCXMS(16). These differences can be attributed to corresponding differences in the pattern of the couplings reported in Table VII. Indeed, the couplings of 2B3u with 1B1g and 1B2u are remarkably larger according to LVCXMS(16). On the contrary, the coupling of 2B3u with 2Ag is larger according to LVCMS(16). The latter also predicts a much larger coupling of the higher-energy state 2B1g with 2Ag explaining why, despite its energy, 2B1g gains some transient population, which, according to LVCMS(16), reaches slightly larger values and decays at a slightly slower rate than in the case of LVCXMS(16).

Analysis of Fig. 4 suggests that after photoexcitation to 1B2u the dynamics is quite simple, being essentially characterized by a progressive (approximatively mono-exponential) flow of population from 1B2u to the lowest-energy state 1B3u. This is not surprising considering that at the FC position, the third state, 1B1g is ∼0.5 eV higher in energy than 1B2u. However, Table VII shows that 1B1g is strongly coupled to both 1B2u and 1B3u states. More specifically, the norm of its coupling to these two states is, respectively, more than three [LVCXMS(16)] and more than four [LVCMS(16)] times larger than the direct 1B1g/1B2u coupling. A small transient population on 1B1g is actually seen in Fig. 4 for the Hamiltonian with the larger couplings [LVCMS(16)]. In Fig. 5, we investigate in greater detail the impact on the 1B2u → 1B3u transfer of the existence of 1B1g and the higher energy states. In order to do that, we compare the dynamics including all the seven coupled states (seven-states model) with a number of reduced models in which some states were removed: the two-state model “1B2u + 1B3u,” the three-state model “1B2u + 1B3u + 1B1g,” and the six-state model obtained including all states except 1B1g. Differences are striking: according to the two-state model, the population transfer is much slower, smaller in amplitude, and shows large oscillations. Including also 1B1g, the population transfer becomes much faster (even more than in the seven-states model) and irreversible, without any significant quantum beating. However, higher energy states also play a role. This is shown considering the six-state model in which 1B1g is removed. In the six-state-model, the predicted population flow from 1B2u to 1B3u is, in fact, similar to what is obtained with the complete seven-state model. Actually, in the long-time limit, 1B3u reaches even a higher population, although the transfer is slower in the first 500 fs (this is better shown by a zoomed-in view of the figure reported in Fig. S10 in the supplementary material). In summary, the existence of 1B1g has a dramatic impact on the 1B2u → 1B3u transfer, much larger than what one could hypothesize looking at the small transient population it acquires. Its main role, in fact, is to provide an alternative and very effective coupling channel between the two lowest states. On the short-time scale, the effect of 1B1g is partially contrasted by the higher-energy states, which slow down the rise of the population of 1B3u. On the long-time scale, however, according to the seven-state model, 1B1g maintains a weak population (∼3%). If such a state is not included in the calculation, this small population flows to 1B3u, making the yield of this state even larger (six-state model). 1B1g and higher-energy states play a qualitatively similar role also according to the XMS(8:16) parameterization, but couplings with 1B1g are smaller. In conclusion, the faster 1B2u → 1B3u decay predicted by LVCMS(16) with respect to LVCXMS(16) is not only due to the larger direct coupling (as discussed above) but also, for a significant part, due to the larger couplings of both states with 1B1g (see Table VII).

FIG. 5.

Dynamics of the populations of the diabatic electronic states after an initial excitation on 1B2u. Comparison of the results obtained with the complete seven-states model and with a number of reduced-dimensionality models in which some electronic states are removed from the LVCMS(16) Hamiltonian.

FIG. 5.

Dynamics of the populations of the diabatic electronic states after an initial excitation on 1B2u. Comparison of the results obtained with the complete seven-states model and with a number of reduced-dimensionality models in which some electronic states are removed from the LVCMS(16) Hamiltonian.

Close modal

Figure 6 plots the diabatic LVC PES at the average position of the WP as a function of time according to the LVCMS(16) Hamiltonian [results for LVCXMS(16) are very similar and are given in Fig. S12 of the supplementary material]. It shows that at all times, S1 and S2 are well separated in energy and rather distant from two pairs of close-lying states, namely, S3–S4, and S5–S6. Interestingly, these data indicate that the average position of the WP does not encounter conical intersections. This finding, together with the smooth changes of the electronic populations, suggests that the picture that better describes the dynamics is not a ballistic movement of the WP toward a CoI. On the contrary, we observe a gradual transfer due to the fact that vibrational states of the upper electronic states are embedded in (and coupled to) a denser manifold of vibrational states of the lower-energy electronic states. Actually, the possible occurrence of fast population transfers in QD even in cases where CoIs are inaccessible has been recently discussed in the literature.63 While this mechanism could be anticipated for an initial excitation to 1B2u, since the initial potential energy of the WP is 3.75 eV (Table II) and the lowest 1B1g/1B2u crossing is at ∼4.2 eV (Table VI), it is noteworthy that the same picture applies also for an initial excitation to 2B3u although several crossings between diabatic states are reachable at this energy, including the (quasi) triple-crossings 1B3u/1B2u/1B1g and 1B1g/2Ag/2B3u.

FIG. 6.

Diabatic LVC potential energies at the average position of the wavepacket obtained for an initial photoexcitation to 1B2u (left) or 2B3u (right) with the LVCMS(16) Hamiltonian. A comparison with the adiabatic energies, very similar, is shown in Fig. S11.

FIG. 6.

Diabatic LVC potential energies at the average position of the wavepacket obtained for an initial photoexcitation to 1B2u (left) or 2B3u (right) with the LVCMS(16) Hamiltonian. A comparison with the adiabatic energies, very similar, is shown in Fig. S11.

Close modal

Finally, Fig. 7 reports the expectation values of all the total-symmetric modes as a function of time for the LVCMS(16) Hamiltonian. The results for LVCXMS(16) are shown in Fig. S13 and are very similar. Both starting from 1B2u and 2B3u, the dynamics is dominated by the oscillations of four modes: two CC stretchings with frequencies 1456 cm−1 (mode 52) and 1669 cm−1 (mode 62) and two lower frequency modes corresponding to a breathing mode with frequency 593 cm−1 (mode 17) and to an in-plane elongation along the long molecular axis with frequency of 406 cm−1 (mode 8). These modes agree with Raman signatures of 1B2u and 2B3u,64,65 and their involvement is consistent with the analysis of excited state vibrational coherences resolved recently in transient absorption spectra with ultrahigh time-resolution (6 fs).39 It is noteworthy that despite the involvement of multiple electronic states coupled differently to the Ag vibrational modes, the dynamics of the average position along individual modes shows only minor deviations in the first 500 fs aside from mode 62, which shows a characteristic shift and damping.

FIG. 7.

Time evolution of average position of the Ag modes for excitation on 1B2u (top) or 2B3u (bottom) excited states. Only the modes with largest displacement are labeled (LVCMS(16) Hamiltonian).

FIG. 7.

Time evolution of average position of the Ag modes for excitation on 1B2u (top) or 2B3u (bottom) excited states. Only the modes with largest displacement are labeled (LVCMS(16) Hamiltonian).

Close modal

We conclude this section mentioning that LVCXMS(12) predicts a very different dynamics (Fig. S14 in the supplementary material), characterized by the fact that both starting from 1B2u and 2B3u, the states 2B3u and 2Ag behave similarly, with very similar populations at all times. Such a peculiar behavior can be explained with the very large coupling between these two states predicted at this level of theory (Table S10 in the supplementary material).

In this contribution, we have combined highly accurate, multiconfigurational electronic structure methods such as RASPT2/RASSCF with a maximum-overlap diabatization technique to parameterize a LVC Hamiltonian for QD. As a case study, we have applied our protocol to the fast QD of pyrene photoexcited to either the first or the second bright state. The rigidity of this molecule justifies the LVC approximation to describe the potential energy surfaces. Yet, its electronic structure and the large number of modes make necessary the inclusion of many electronic states and the development of an effective diabatization protocol to build the vibronic Hamiltonian. From the point of view of the electronic structure theory, several characteristics of pyrene require the adoption of multiconfigurational methods, such as the presence of a state with a high contribution from a double excitation (2Ag) and the difficulty of many TD-DFT functionals in reproducing the relative order of the lowest-energy states.66–68 The involvement in the dynamics of the 2Ag state makes also problematic the usage of methods such as ADC(2) or CC2 since although they have shown remarkable accuracy for single excitations in organic molecules, they are not accurate for double-excited states.69,70 The parameterization based on RASPT2/RASSCF makes also our LVC Hamiltonian suitable for the simulation, in the near future, of transient absorption spectra. To this end, in fact, the computation of transition dipoles with the possible final states reached by the absorption of the probe is required, which have an increased probability to show a significant double-excited character.71 To the best of our knowledge, this is the first reported example of LVC parameterization based on energies and WFs overlaps computed with RASPT2/RASSCF electronic structure calculations. Our results evidence that it is not a “black box” procedure. While, in principle, the RASTP2/RASSCF protocol is able to describe states with different nature on an equal footing, large active spaces, beyond the full-π set of orbitals, are needed to achieve this. Therefore, benchmarking is essential for assuring the convergence of the excited state energies with respect to the active space size.72 The undertaking is nowadays possible even for relatively big systems, thanks to flexible approaches to the construction of the active space such as the generalized active space (GAS)SCF/GASPT2 approach73 or the generalized multi-configuration quasi-degenerate perturbation theory (GMCQDPT),74 as well as modern day CI solvers such as the density matrix renormalization group (DMRG)75 and the full configuration interaction quantum Monte Carlo (FCIQMC),76 to name a few, which allow to handle active spaces with many tens of orbitals. Another critical point to address is the various flavors of the perturbative correction each one with its strengths and weaknesses. Our result indicates that SS-RASPT2 should be used for the energy calculations whenever the electronic states are far apart in energy. On the other hand, MS- and XMS-RASPT2 energies and WFs more reliable in the case of close-lying, interacting states. In particular, perturbatively modified WFs should be used in the maximum-overlap diabatization procedure. Finally, when symmetry can be applied to reduce the computational cost, attention is advised in regard to biases introduced by the RASSCF/RASPT2 protocol in the calculation of coupling parameters along symmetry-reducing normal modes. Exemplarily, the unbalanced description of two close lying states (i.e., 2Ag, already well described with the full-π active space, and the 2B3u that shows a strong dependence on the active space size) could result in non-physically large vibronic couplings as demonstrated by LVC parameterization at the XMS(8:12) level.

In spite of these complications, benchmarking of diabatic PESs obtained with our model [parameterized at adequate levels such as MS(8:16) and XMS(8:16)] against RASPT2 calculations proves that the LVC Hamiltonian can be highly accurate, being also able to predict the structure and energy of both excited state minima and crossings between the states included in the model. LVCXMS(16) and LVCMS(16) dynamics are qualitatively similar. Still, both for an initial excitation to 1B2u and to 2B3u, LVCMS(16) predicts that the decay from 1B2u to 1B3u is remarkably faster. These differences point out that, at the state of the art, even quite sophisticated electronic structure methods cannot guarantee the computation of precise decay times. On the one side, this result witnesses the necessity to use accurate methods even for the parameterization of simple vibronic Hamiltonians such as LVC. On the other side, it documents the necessity of further efforts in the development of electronic structure methods for excited states of medium size molecules.

The QD simulations indicate that after an initial photoexcitation to 1B2u (S2), the population progressively flows to 1B3u (S1). In particular, the population growth with a sub-100 fs time constant predicted by the LVCMS(16) Hamiltonian agrees very well with experimental observations.34,36 Quite interestingly, this transfer is strongly affected by the existence of higher-energy states, especially 1B1g, even if it lies ∼0.5 eV above the bright state in the Franck–Condon region. This finding highlights that in order to obtain robust QD results, it is necessary to adopt LVC models including a sufficiently large number of diabatic states. Direct excitation of the second bright state 2B3u (S5) leads to its ultrafast (sub-100 fs) depopulation in favor of a number of intermediate states, especially 1B2u, followed by a much slower progressively decay to 1B3u (S1), supporting the mechanism proposed based on recent experimental findings.36 Rather surprisingly, in both QD simulations, population transfers occur smoothly and in an ultrafast manner even if the average position of the WP never get really close to crossing points of the diabatic (and adiabatic) states. In particular, the 1B2u (S2) → 1B3u (S1) transfer was found to occur on a sub-100 fs time-scale even if the CoI lies ∼0.4 eV above the FC point. This observation can be rationalized by coupling between vibrational levels, rather than ballistic motion toward a CoI. In the light of this finding, the question arises whether semi-classical trajectory-based approaches, which treat nuclei classically, are capable of capturing the ultrafast nature of the internal conversion.

Finally, it is noteworthy that the protocol for the parametrization of LVC Hamiltonians from RASPT2/RASSCF is fully general and ready to be applied to other interesting problems, like the ultrafast internal conversion in photoexcited nucleobases.19 Furthermore, the protocol is straight-forwardly extendable to incorporate spin–orbit couplings to describe inter-system crossing.77 

See the supplementary material for pyrene normal modes and frequencies, adiabatic excited state minima with the LVC displaced harmonic oscillator model, adiabatic overlap matrices, ML-MCTDH trees and convergence tests, population dynamics of models with reduced number of electronic states, diabatic and adiabatic energies for the diabatic states minima and conical intersections estimated by LVC, diabatic and adiabatic potential energy surfaces at the average position of the wavepacket for the LVCMS(16) and LVCXMS(16) Hamiltonians, average position of the Ag modes during wavepacket propagations for the LVCXMS(16) Hamiltonian, comparison of the population dynamics with the LVCMS(16), LVCXMS(12), and LVCXMS(16) Hamiltonians, and norm of the diabatic coupling vectors for LVCXMS(12) parameterization.

F.A. and D.A. contributed equally to this work.

This project has received funding from the European Union’s Horizon 2020 Research and Innovation Programme under Marie Sklodowska-Curie Grant Agreement No. 765266 (LightDyNAmics). D.A. acknowledges Fundación Ramón Areces (Spain) for funding his postdoctoral stay at ICCOM-CNR Pisa.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
C.
Leforestier
,
R. H.
Bisseling
,
C.
Cerjan
,
M. D.
Feit
,
R.
Friesner
,
A.
Guldberg
,
A.
Hammerich
,
G.
Jolicard
,
W.
Karrlein
,
H.-D.
Meyer
,
N.
Lipkin
,
O.
Roncero
, and
R.
Kosloff
, “
A comparison of different propagation schemes for the time dependent Schrödinger equation
,”
J. Comput. Phys.
94
,
59
80
(
1991
).
2.
M. H.
Beck
,
A.
Jäckle
,
G. A.
Worth
, and
H.-D.
Meyer
, “
The multiconfiguration time-dependent Hartree (MCTDH) method: A highly efficient algorithm for propagating wavepackets
,”
Phys. Rep.
324
,
1
105
(
2000
).
3.
Multidimensional Quantum Dynamics: MCTDH Theory and Applications
, edited by
H.-D.
Meyer
,
F.
Gatti
, and
G. A.
Worth
(
WileyVCH
,
Weinheim
,
2009
).
4.
F.
Gatti
,
B.
Lasorne
,
H.-D.
Meyer
, and
A.
Nauts
,
Applications of Quantum Dynamics in Chemistry
, Lecture Notes in Chemistry, Vol. 98 (
Springer International Publishing
,
2017
).
5.
H.
Wang
and
M.
Thoss
, “
Multilayer formulation of the multiconfiguration time-dependent Hartree theory
,”
J. Chem. Phys.
119
,
1289
1299
(
2003
).
6.
U.
Manthe
, “
A multilayer multiconfigurational time-dependent Hartree approach for quantum dynamics on general potential energy surfaces
,”
J. Chem. Phys.
128
,
164116
(
2008
).
7.
U.
Manthe
, “
Layered discrete variable representations and their application within the multiconfigurational time-dependent Hartree approach
,”
J. Chem. Phys.
130
,
054109
(
2009
).
8.
O.
Vendrell
and
H.-D.
Meyer
, “
Multilayer multiconfiguration time-dependent Hartree method: Implementation and applications to a Henon–Heiles Hamiltonian and to pyrazine
,”
J. Chem. Phys.
134
,
044135
(
2011
).
9.
H.
Köuppel
,
W.
Domcke
, and
L. S.
Cederbaum
, “
Multimode molecular dynamics beyond the Born-Oppenheimer approximation
,”
Adv. Chem. Phys.
57
,
59
(
1984
).
10.
H.
Köppel
, “
Diabatic representation: Methods for the construction of diabatic electronic state
,” in
Conical Intersections, Electronic Structure, Dynamics and Spectroscopy
, edited by
W.
Domcke
,
R.
Yarkony
, and
H.
Köppel
(
World Scientific Publishing Co.
,
Singapore
,
2004
), pp.
175
204
.
11.
G.
Herzberg
and
H. C.
Longuet-Higgins
, “
Intersection of potential energy surfaces in polyatomic molecules
,”
Discuss. Faraday Soc.
35
,
77
82
(
1963
).
12.
T. J.
Penfold
,
E.
Gindensperger
,
C.
Daniel
, and
C. M.
Marian
, “
Spin-vibronic mechanism for intersystem crossing
,”
Chem. Rev.
118
(
15
),
6975
7025
(
2018
).
13.
M.
Fumanal
,
E.
Gindensperger
, and
C.
Daniel
, “
Ultrafast intersystem crossing vs internal conversion in α-diimine transition metal complexes: Quantum evidence
,”
J. Phys. Chem. Lett.
9
(
17
),
5189
5195
(
2018
).
14.
F.
Plasser
,
S.
Mai
,
M.
Fumanal
,
E.
Gindensperger
,
C.
Daniel
, and
L.
González
, “
Strong Influence of Decoherence Corrections and Momentum Rescaling in Surface Hopping Dynamics of Transition Metal Complexes
,”
J. Chem. Theory Comput.
15
(
9
),
5031
5045
(
2019
).
15.
S.
Mai
and
L.
González
, “
Identification of important normal modes in nonadiabatic dynamics simulations by coherence, correlation, and frequency analyses
,”
J. Chem. Phys.
151
(
24
),
244115
(
2019
).
16.
D.
Picconi
,
A.
Lami
, and
F.
Santoro
, “
Hierarchical transformation of Hamiltonians with linear and quadratic couplings for nonadiabatic quantum dynamics: Application to the ππ*/nπ* internal conversion in thymine
,”
J. Chem. Phys.
136
,
244104
(
2012
).
17.
D.
Picconi
,
F. J.
Avila Ferrer
,
R.
Improta
,
A.
Lami
, and
F.
Santoro
, “
Quantum-classical effective-modes dynamics of the ππ* → nπ* decay in 9H-adenine. A quadratic vibronic coupling model
,”
Faraday Discuss.
163
,
223
242
(
2013
).
18.
J.
Cerezo
,
Y.
Liu
,
N.
Lin
,
X.
Zhao
,
R.
Improta
, and
F.
Santoro
, “
Mixed quantum/classical method for nonadiabatic quantum dynamics in explicit solvent models: The ππ*/nπ* decay of thymine in water as a test case
,”
J. Chem. Theory Comput.
14
,
820
832
(
2018
).
19.
M.
Yaghoubi Jouybari
,
Y.
Liu
,
R.
Improta
, and
F.
Santoro
, “
The ultrafast dynamics of the two lowest bright excited states of cytosine and 1-methyl-cytosine: A quantum dynamical study
,”
J. Chem. Theory Comput.
16
(
9
),
5792
5808
(
2020
).
20.
M.
Yaghoubi Jouybari
,
Y.
Liu
,
R.
Improta
, and
F.
Santoro
, “
Quantum dynamics of the ππ*/nπ* decay of the epigenetic nucleobase 1,5-dimethyl-cytosine in the gas phase
,”
Phys. Chem. Chem. Phys.
22
,
26525
26535
(
2020
).
21.
M.
Pápai
,
M.
Abedi
,
G.
Levi
,
E.
Biasin
,
M. M.
Nielsen
, and
K. B.
Møller
, “
Theoretical evidence of solvent-mediated excited-state dynamics in a functionalized iron sensitizer
,”
J. Phys. Chem. C
123
(
4
),
2056
2065
(
2019
).
22.
R.
Improta
,
V.
Barone
,
A.
Lami
, and
F.
Santoro
, “
Quantum dynamics of the ultrafast ππ*/nπ* population transfer in uracil and 5-fluoro-uracil in water and acetonitrile
,”
J. Phys. Chem. B
113
,
14491
14503
(
2009
).
23.
J.
Seibt
,
T.
Winkler
,
K.
Renziehausen
,
V.
Dehm
,
F.
Würthner
,
H.-D.
Meyer
, and
V.
Engel
, “
Vibronic transitions and quantum dynamics in molecular oligomers: A theoretical analysis with an application to aggregates of perylene bisimides
,”
J. Phys. Chem. A
113
(
48
),
13475
13482
(
2009
).
24.
H.
Tamura
,
I.
Burghardt
, and
M.
Tsukada
, “
Exciton dissociation at thiophene/fullerene interfaces: The electronic structures and quantum dynamics
,”
J. Phys. Chem. C
115
(
20
),
10205
10210
(
2011
).
25.
M.
Keβ
,
G.
Worth
, and
V.
Engel
, “
Two-dimensional vibronic spectroscopy of molecular aggregates: Trimers, dimers, and monomers
,”
J. Chem. Phys.
145
(
8
),
084305
(
2016
).
26.
F.
Segatta
,
L.
Cupellini
,
M.
Garavelli
, and
B.
Mennucci
, “
Quantum chemical modeling of the photoinduced activity of multichromophoric biosystems
,”
Chem. Rev.
119
(
16
),
9361
9380
(
2019
).
27.
D.
Aranda
and
F.
Santoro
, “
Vibronic spectra of π-conjugated systems with a multitude of coupled states. A protocol based on linear vibronic coupling models and quantum dynamics tested on hexahelicene
,”
J. Chem. Theory Comput.
(published online).
28.
R.
Cimiraglia
,
J.-P.
Malrieu
,
M.
Persico
, and
F.
Spiegelmann
, “
Quasi-diabatic states and dynamical couplings from ab initio CI calculations: A new proposal
,”
J. Phys. B: At. Mol. Phys.
18
(
15
),
3073
(
1985
).
29.
J.
Finley
,
P.-Å.
Malmqvist
,
B. O.
Roos
, and
L.
Serrano-Andrés
, “
The multi-state CASPT2 method
,”
Chem. Phys. Lett.
288
(
2-4
),
299
306
(
1998
).
30.
T.
Shiozaki
,
W.
Győrffy
,
P.
Celani
, and
H.-J.
Werner
, “
Communication: Extended multi-state complete active space second-order perturbation theory: Energy and nuclear gradients
,”
J. Chem. Phys.
135
(
8
),
081106
(
2011
).
31.
S.
Battaglia
and
R.
Lindh
, “
Extended dynamically weighted CASPT2: The best of two worlds
,”
J. Chem. Theory Comput.
16
(
3
),
1555
1567
(
2020
).
32.
J. W.
Park
, “
Single-state single-reference and multistate multireference zeroth-order Hamiltonians in MS-CASPT2 and conical intersections
,”
J. Chem. Theory Comput.
15
(
7
),
3960
3973
(
2019
).
33.
P.
Foggi
,
L.
Pettini
,
I.
Santa
,
R.
Righini
, and
S.
Califano
, “
Transient absorption and vibrational relaxation dynamics of the lowest excited singlet state of pyrene in solution
,”
J. Phys. Chem.
99
(
19
),
7439
7445
(
1995
).
34.
F. V. R.
Neuwahl
and
P.
Foggi
, “
Direct observation of s2–s1 internal conversion in pyrene by femtosecond transient absorption
,”
Laser Chem.
19
(
1-4
),
375
379
(
1999
).
35.
M.
Raytchev
,
E.
Pandurski
,
I.
Buchvarov
,
C.
Modrakowski
, and
T.
Fiebig
, “
Bichromophoric interactions and time-dependent excited state mixing in pyrene derivatives. A femtosecond broad-band pump-probe study
,”
J. Phys. Chem. A
107
(
23
),
4592
4600
(
2003
).
36.
R.
Borrego-Varillas
,
L.
Ganzer
,
G.
Cerullo
, and
C.
Manzoni
, “
Ultraviolet transient absorption spectrometer with sub-20-fs time resolution
,”
Appl. Sci.
8
(
6
),
989
(
2018
).
37.
A.
Nenov
,
A.
Giussani
,
B. P.
Fingerhut
,
I.
Rivalta
,
E.
Dumont
,
S.
Mukamel
, and
M.
Garavelli
, “
Spectral lineshapes in nonlinear electronic spectroscopy
,”
Phys. Chem. Chem. Phys.
17
(
46
),
30925
30936
(
2015
).
38.
M. K.
Roos
,
S.
Reiter
, and
R.
de Vivie-Riedle
, “
Ultrafast relaxation from 1La to 1Lb in pyrene: A theoretical study
,”
Chem. Phys.
515
,
586
595
(
2018
).
39.
A.
Picchiotti
,
A.
Nenov
,
A.
Giussani
,
V. I.
Prokhorenko
,
R. J.
Dwayne Miller
,
S.
Mukamel
, and
M.
Garavelli
, “
Pyrene, a test case for deep-ultraviolet molecular photophysics
,”
J. Phys. Chem. Lett.
10
(
12
),
3481
3487
(
2019
).
40.
J. A.
Noble
,
C.
Aupetit
,
D.
Descamps
,
S.
Petit
,
A.
Simon
,
J.
Mascetti
,
N.
Ben Amor
, and
V.
Blanchet
, “
Ultrafast electronic relaxations from the S3 state of pyrene
,”
Phys. Chem. Chem. Phys.
21
(
26
),
14111
14125
(
2019
).
41.
J.
Neugebauer
,
E. J.
Baerends
, and
M.
Nooijen
, “
Vibronic coupling and double excitations in linear response time-dependent density functional calculations: Dipole-allowed states of N2
,”
J. Chem. Phys.
121
(
13
),
6155
6166
(
2004
).
42.
Y.
Liu
,
L.
Martínez-Fernández
,
J.
Cerezo
,
G.
Prampolini
,
R.
Improta
, and
F.
Santoro
, “
Multistate coupled quantum dynamics of photoexcited cytosine in gas-phase: Nonadiabatic absorption spectrum and ultrafast internal conversions
,”
Chem. Phys.
515
,
452
463
(
2018
).
43.
H.
Köppel
,
W.
Domcke
, and
L.
Cederbaum
, “
The multi-mode vibronic-coupling approach
,” in
Conical Intersections, Electronic Structure, Dynamics and Spectroscopy
, edited by
W.
Domcke
,
R.
Yarkony
, and
H.
Köppel
(
World Scientific Publishing Co.
,
Singapore
,
2004
), pp.
323
368
.
44.
F.
Plasser
,
S.
Gómez
,
M. F. S. J.
Menger
,
S.
Mai
, and
L.
González
, “
Highly efficient surface hopping dynamics using a linear vibronic coupling model
,”
Phys. Chem. Chem. Phys.
21
,
57
69
(
2019
).
45.
B.
Gonon
,
B.
Lasorne
,
G.
Karras
,
L.
Joubert-Doriol
,
D.
Lauvergnat
,
F.
Billard
,
B.
Lavorel
,
O.
Faucher
,
S.
Guérin
,
E.
Hertz
, and
F.
Gatti
, “
A generalized vibronic-coupling Hamiltonian for molecules without symmetry: Application to the photoisomerization of benzopyran
,”
J. Chem. Phys.
150
(
12
),
124109
(
2019
).
46.
M.
Sala
,
B.
Lasorne
,
F.
Gatti
, and
S.
Guérin
, “
The role of the low-lying dark nπ* states in the photophysics of pyrazine: A quantum dynamics study
,”
Phys. Chem. Chem. Phys.
16
,
15957
15967
(
2014
).
47.
A.
Lehr
,
S.
Gómez
,
M. A.
Parkes
, and
G. A.
Worth
, “
The role of vibronic coupling in the electronic spectroscopy of maleimide: A multi-mode and multi-state quantum dynamics study
,”
Phys. Chem. Chem. Phys.
22
,
25272
25283
(
2020
).
48.
V.
Sauri
,
L.
Serrano-Andrés
,
A. R. M.
Shahi
,
L.
Gagliardi
,
S.
Vancoillie
, and
K.
Pierloot
, “
Multiconfigurational second-order perturbation theory restricted active space (RASPT2) method for electronic excited states: A benchmark study
,”
J. Chem. Theory Comput.
7
(
1
),
153
168
(
2010
).
49.
A.
Nenov
,
S.
Mukamel
,
M.
Garavelli
, and
I.
Rivalta
, “
Two-dimensional electronic spectroscopy of benzene, phenol, and their dimer: An efficient first-principles simulation protocol
,”
J. Chem. Theory Comput.
11
(
8
),
3755
3771
(
2015
).
50.
A.
Nenov
,
A.
Giussani
,
J.
Segarra-Martí
,
V. K.
Jaiswal
,
I.
Rivalta
,
G.
Cerullo
,
S.
Mukamel
, and
M.
Garavelli
, “
Modeling the high-energy electronic state manifold of adenine: Calibration for nonlinear electronic spectroscopy
,”
J. Chem. Phys.
142
(
21
),
212443
(
2015
).
51.
A.
Nenov
,
I.
Rivalta
,
S.
Mukamel
, and
M.
Garavelli
, “
Bidimensional electronic spectroscopy on indole in gas phase and in water from first principles
,”
Comput. Theor. Chem.
1040-1041
,
295
303
(
2014
).
52.
I.
Fdez. Galván
,
M.
Vacher
,
A.
Alavi
,
C.
Angeli
,
F.
Aquilante
,
J.
Autschbach
,
J. J.
Bao
,
S. I.
Bokarev
,
N. A.
Bogdanov
,
R. K.
Carlson
,
L. F.
Chibotaru
,
J.
Creutzberg
,
N.
Dattani
,
M. G.
Delcey
,
S.
Dong
,
A.
Dreuw
,
L.
Freitag
,
L.
Manuel Frutos
,
L.
Gagliardi
,
F.
Gendron
,
A.
Giussani
,
L.
González
,
G.
Grell
,
M.
Guo
,
C. E.
Hoyer
,
M.
Johansson
,
S.
Keller
,
S.
Knecht
,
G.
Kovačević
,
E.
Källman
,
G.
Li Manni
,
M.
Lundberg
,
Y.
Ma
,
S.
Mai
,
J.
Pedro Malhado
,
P. Å.
Malmqvist
,
P.
Marquetand
,
S. A.
Mewes
,
J.
Norell
,
M.
Olivucci
,
M.
Oppel
,
Q. M.
Phung
,
K.
Pierloot
,
F.
Plasser
,
M.
Reiher
,
A. M.
Sand
,
I.
Schapiro
,
P.
Sharma
,
C. J.
Stein
,
L. K.
Sørensen
,
D. G.
Truhlar
,
M.
Ugandi
,
L.
Ungur
,
A.
Valentini
,
S.
Vancoillie
,
V.
Veryazov
,
O.
Weser
,
T. A.
Wesołowski
,
P.-O.
Widmark
,
S.
Wouters
,
A.
Zech
,
J. P.
Zobel
, and
R.
Lindh
, “
OpenMolcas: From source code to insight
,”
J. Chem. Theory Comput.
15
(
11
),
5925
5964
(
2019
).
53.
F.
Aquilante
,
J.
Autschbach
,
A.
Baiardi
,
S.
Battaglia
,
V. A.
Borin
,
L. F.
Chibotaru
,
I.
Conti
,
L.
De Vico
,
M.
Delcey
,
I.
Fdez. Galván
,
N.
Ferré
,
L.
Freitag
,
M.
Garavelli
,
X.
Gong
,
S.
Knecht
,
E. D.
Larsson
,
R.
Lindh
,
M.
Lundberg
,
P. Å.
Malmqvist
,
A.
Nenov
,
J.
Norell
,
M.
Odelius
,
M.
Olivucci
,
T. B.
Pedersen
,
L.
Pedraza-González
,
Q. M.
Phung
,
K.
Pierloot
,
M.
Reiher
,
I.
Schapiro
,
J.
Segarra-Martí
,
F.
Segatta
,
L.
Seijo
,
S.
Sen
,
D.-C.
Sergentu
,
C. J.
Stein
,
L.
Ungur
,
M.
Vacher
,
A.
Valentini
, and
V.
Veryazov
, “
Modern quantum chemistry with [OPEN]Molcas
,”
J. Chem. Phys.
152
(
21
),
214117
(
2020
).
54.
G. A.
Worth
,
K.
Giri
,
G. W.
Richings
,
M. H.
Beck
,
A.
Jäckle
, and
H.-D.
Meyer
, The QUANTICS Package, version 1.1,
University of Birmingham
,
Birmingham, UK
(
2015
).
55.
G. A.
Worth
, “
QUANTICS: A general purpose package for quantum molecular dynamics simulations
,”
Comput. Phys. Commun.
248
,
107040
(
2020
).
56.
J.
Ferguson
,
L. W.
Reeves
, and
W. G.
Schneider
, “
Vapor absorption spectra of naphthalene, anthracene and pyrene
,”
Can. J. Chem.
35
(
10
),
1117
1136
(
1957
).
57.
H.
Baba
and
M.
Aoi
, “
Vapor-phase fluorescence spectra from the second excited singlet state of pyrene and its derivatives
,”
J. Mol. Spectrosc.
46
(
2
),
214
222
(
1973
).
58.
P. R.
Salvi
,
P.
Foggi
, and
E.
Castellucci
, “
The two-photon excitation spectrum of pyrene
,”
Chem. Phys. Lett.
98
(
3
),
206
211
(
1983
).
59.
D. S.
Karpovich
and
G. J.
Blanchard
, “
Relating the polarity-dependent fluorescence response of pyrene to vibronic coupling. Achieving a fundamental understanding of the py polarity scale
,”
J. Phys. Chem.
99
(
12
),
3951
3958
(
1995
).
60.
A.
Ya. Freidzon
,
R. R.
Valiev
, and
A. A.
Berezhnoy
, “
Ab initio simulation of pyrene spectra in water matrices
,”
RSC Adv.
4
(
79
),
42054
42065
(
2014
).
61.
S.
Sen
and
I.
Schapiro
, “
A comprehensive benchmark of the XMS-CASPT2 method for the photochemistry of a retinal chromophore model
,”
Mol. Phys.
116
(
19-20
),
2571
2582
(
2018
).
62.
In the cases were the RASPT2 minimum was not reported in the literature or had been obtained with a different state-averaging (i.e., for states S3, S4, and S6), we have done the SS(8:8) optimization.
63.
C. A.
Farfan
and
D. B.
Turner
, “
Nonadiabatic photochemistry induced by inaccessible conical intersections
,”
J. Phys. Chem. A
123
(
36
),
7768
7776
(
2019
).
64.
C. M.
Jones
and
S. A.
Asher
, “
Ultraviolet resonance Raman study of the pyrene S4, S3, and S2 excited electronic states
,”
J. Chem. Phys.
89
(
5
),
2649
2661
(
1988
).
65.
H.
Shinohara
,
Y.
Yamakita
, and
K.
Ohno
, “
Raman spectra of polycyclic aromatic hydrocarbons. Comparison of calculated Raman intensity distributions with observed spectra for naphthalene, anthracene, pyrene, and perylene
,”
J. Mol. Struct.
442
(
1-3
),
221
234
(
1998
).
66.
E. L.
Graef
and
J. B. L.
Martins
, “
Analysis of lowest energy transitions at TD-DFT of pyrene in vacuum and solvent
,”
J. Mol. Model.
25
(
7
),
183
(
2019
).
67.
B.
Shi
,
D.
Nachtigallová
,
A. J. A.
Aquino
,
F. B. C.
Machado
, and
H.
Lischka
, “
High-level theoretical benchmark investigations of the UV-vis absorption spectra of paradigmatic polycyclic aromatic hydrocarbons as models for graphene quantum dots
,”
J. Chem. Phys.
150
(
12
),
124302
(
2019
).
68.
F. J.
Avila Ferrer
,
V.
Barone
,
C.
Cappelli
, and
F.
Santoro
, “
Duschinsky, Herzberg–Teller, and multiple electronic resonance interferential effects in resonance Raman spectra and excitation profiles. The case of pyrene
,”
J. Chem. Theory Comput.
9
(
8
),
3597
3611
(
2013
).
69.
D.
Jacquemin
,
I.
Duchemin
, and
X.
Blase
, “
0–0 energies using hybrid schemes: Benchmarks of TD-DFT, CIS(D), ADC(2), CC2, and BSE/GW formalisms for 80 real-life compounds
,”
J. Chem. Theory Comput.
11
(
11
),
5340
5359
(
2015
).
70.
P.-F.
Loos
,
F.
Lipparini
,
M.
Boggio-Pasqua
,
A.
Scemama
, and
D.
Jacquemin
, “
A mountaineering strategy to excited states: Highly accurate energies and benchmarks for medium sized molecules
,”
J. Chem. Theory Comput.
16
(
3
),
1711
1741
(
2020
).
71.
J.
Segarra-Martí
,
S.
Mukamel
,
M.
Garavelli
,
A.
Nenov
, and
I.
Rivalta
, “
Towards accurate simulation of two-dimensional electronic spectroscopy
,” in
Topics in Current Chemistry Collections
(
Springer International Publishing
,
2018
), pp.
63
112
.
72.
S.
Shirai
and
S.
Inagaki
, “
Ab initio study on the excited states of pyrene and its derivatives using multi-reference perturbation theory methods
,”
RSC Adv.
10
(
22
),
12988
12998
(
2020
).
73.
D.
Ma
,
G.
Li Manni
,
J.
Olsen
, and
L.
Gagliardi
, “
Second-order perturbation theory for generalized active space self-consistent-field wave functions
,”
J. Chem. Theory Comput.
12
(
7
),
3208
3213
(
2016
).
74.
H.
Nakano
,
R.
Uchiyama
, and
K.
Hirao
, “
Quasi-degenerate perturbation theory with general multiconfiguration self-consistent field reference functions
,”
J. Comput. Chem.
23
(
12
),
1166
1175
(
2002
).
75.
L.
Freitag
and
M.
Reiher
,
The Density Matrix Renormalization Group for Strong Correlation in Ground and Excited States
(
Wiley
,
2020
).
76.
R. J.
Anderson
,
T.
Shiozaki
, and
G. H.
Booth
, “
Efficient and stochastic multireference perturbation theory for large active spaces within a full configuration interaction quantum Monte Carlo framework
,”
J. Chem. Phys.
152
(
5
),
054101
(
2020
).
77.
J.
Eng
,
C.
Gourlaouen
,
E.
Gindensperger
, and
C.
Daniel
, “
Spin-vibronic quantum dynamics for ultrafast excited-state processes
,”
Acc. Chem. Res.
48
(
3
),
809
817
(
2015
).

Supplementary Material