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 (S_{2} or S_{5}). 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 S_{5}, whose ultrafast deactivation mechanism was not yet clearly understood.

## I. INTRODUCTION

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) functionals^{19,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 weighted^{31} 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 experimentally^{33–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.

## II. METHODOLOGY: THE LINEAR VIBRONIC COUPLING MODEL

We consider a *n* dimensional diabatic basis, |**d**⟩ = (|*d*_{1}⟩, |*d*_{2}⟩, …, |*d*_{n}⟩), and the following expression of the Hamiltonian:

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:

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} (*i* ≠ *j*).

The LVC Hamiltonian is parameterized by defining diabatic states |*d*_{i}⟩ to be coincident with the adiabatic reference states |*a*_{i}⟩ 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**)⟩,

The transformation matrix **D** that defines the diabatic states at Δ_{α},

is then obtained as

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(\Delta \alpha )=diag(E1ad(\Delta \alpha ),E2ad(\Delta \alpha ),\u2026,Enad(\Delta \alpha ))$ and the diabatic potential terms are simply

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

In the following, the normal coordinates **q** and frequencies **Ω** were obtained at the second order perturbation theory level (MP2), whereas the energies $Eiad(\Delta \alpha )$ 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}|*d*_{i}⟩|Ψ_{i}(**q**, *t*)⟩, and the time evolution is computed by solving the time-dependent Schrödinger equation,

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 *P*_{i}(*t*) = ⟨Ψ_{i}(**q**, *t*)|Ψ_{i}(**q**, *t*)⟩.

## III. COMPUTATIONAL DETAILS

### A. Electronic structure calculations

Pyrene is a highly symmetric molecule (D_{2h} 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 S_{0} equilibrium geometry, belonging to four different irreducible representations: one state in A_{g}, two in B_{3u}, one in B_{2u}, and three in B_{1g}. 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 $\u27e8Siref|Sjdispl\u27e9$ 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 (A_{g}) modes, interstate couplings also exist along modes belonging to B_{1g}, B_{2u}, and B_{3u} 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 D_{2h} symmetry labels of the adiabatic states they coincide with at the S_{0} minimum. Adiabatic states, on the contrary, will be denoted with the usual nomenclature S_{x} 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.

Irreducible representation . | Point group at . | Classification of D_{2h} states
. | |
---|---|---|---|

of modes . | displaced geometries . | into new irreducible representations . | |

B_{3u} | C_{2v} | A_{1}: 1A_{g}, 1B_{3u}, 2A_{g}, 2B_{3u} | B_{1}: 1B_{2u}, 1B_{1g}, 2B_{1g}, 3B_{1g} |

B_{2u} | C_{2v} | A_{1}: 1A_{g}, 1B_{2u}, 2A_{g} | B_{2}: 1B_{3u}, 1B_{1g}, 2B_{3u}, 2B_{1g}, 3B_{1g} |

B_{1g} | C_{2h} | A_{g}: 1A_{g}, 2A_{g}, 1B_{1g}, 2B_{1g}, 3B_{1g} | B_{u}: 1B_{3u}, 1B_{2u}, 2B_{3u} |

Irreducible representation . | Point group at . | Classification of D_{2h} states
. | |
---|---|---|---|

of modes . | displaced geometries . | into new irreducible representations . | |

B_{3u} | C_{2v} | A_{1}: 1A_{g}, 1B_{3u}, 2A_{g}, 2B_{3u} | B_{1}: 1B_{2u}, 1B_{1g}, 2B_{1g}, 3B_{1g} |

B_{2u} | C_{2v} | A_{1}: 1A_{g}, 1B_{2u}, 2A_{g} | B_{2}: 1B_{3u}, 1B_{1g}, 2B_{3u}, 2B_{1g}, 3B_{1g} |

B_{1g} | C_{2h} | A_{g}: 1A_{g}, 2A_{g}, 1B_{1g}, 2B_{1g}, 3B_{1g} | B_{u}: 1B_{3u}, 1B_{2u}, 2B_{3u} |

Electronic structure calculations with D_{2h} 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, 8$\u2225$0, 0$\u2225$4, 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, 8$\u2225$0, 0$\u2225$4, 12) and RAS(4, 8$\u2225$0, 0$\u2225$4, 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–51} Figure 1 shows the active orbitals.

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 D_{2h} symmetry (at the reference geometry and along A_{g} 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 B_{1g} 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 D_{2h} 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 A_{1} states at geometries with C_{2v} symmetry obtained by displacing along B_{3u} 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.

### B. QD calculations

ML-MCTDH wavepacket propagation^{2–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.

## IV. RESULTS AND DISCUSSION

### A. Energy calculations

The lowest seven excited states of pyrene belong to four irreducible representations (Table II). Among these states, we identify two optically bright states—1B_{2u} with dominant configuration H(OMO) → L(UMO) and 2B_{3u} 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 A_{g} symmetry in the vicinity of the second bright state evidencing the need of multiconfigurational methods.

. | . | . | TDM . | Energy (eV) . | Experimental . | ||||
---|---|---|---|---|---|---|---|---|---|

State . | Label . | CSFs . | (Debye) . | SS(8:8) . | SS(8:12) . | SS(8:16) . | MS(8:16) . | XMS(8:16) . | ΔE_{0-0} (eV)
. |

S_{0} | 1A_{g} | GS | … | 0.00 | 0.00 | 0.00 | … | … | … |

S_{1} | 1B_{3u} | H → L + 1 | 0.00 | 3.23 | 3.22 | 3.23 | … | … | 3.36^{57} |

H-1 → L | |||||||||

S_{2} | 1B_{2u} | H → L | 1.83 | 3.55 | 3.69 | 3.75 | … | … | 3.84^{56} |

S_{3} | 1B_{1g} | H → L + 2 | 0.00 | 4.11 | 4.13 | 4.16 | 4.00 | 4.10 | 4.12^{58} |

S_{4} | 2A_{g} | (H → L)^{2} | 0.00 | 4.30 | 4.35 | 4.32 | … | … | 4.29^{58} |

S_{5} | 2B_{3u} | H → L + 1 | 1.73 | 4.18 | 4.35 | 4.43 | … | … | 4.66^{56} |

H-1 → L | |||||||||

S_{6} | 2B_{1g} | H-2 → L | 0.00 | 4.28 | 4.46 | 4.56 | 4.64 | 4.48 | 4.54^{58} |

H → L + 2 | |||||||||

S_{7} | 3B_{1g} | H-3 → L | 0.00 | 4.73 | 4.77 | 4.82 | 4.89 | 4.85 | 4.94^{58} |

H → L + 3 |

. | . | . | TDM . | Energy (eV) . | Experimental . | ||||
---|---|---|---|---|---|---|---|---|---|

State . | Label . | CSFs . | (Debye) . | SS(8:8) . | SS(8:12) . | SS(8:16) . | MS(8:16) . | XMS(8:16) . | ΔE_{0-0} (eV)
. |

S_{0} | 1A_{g} | GS | … | 0.00 | 0.00 | 0.00 | … | … | … |

S_{1} | 1B_{3u} | H → L + 1 | 0.00 | 3.23 | 3.22 | 3.23 | … | … | 3.36^{57} |

H-1 → L | |||||||||

S_{2} | 1B_{2u} | H → L | 1.83 | 3.55 | 3.69 | 3.75 | … | … | 3.84^{56} |

S_{3} | 1B_{1g} | H → L + 2 | 0.00 | 4.11 | 4.13 | 4.16 | 4.00 | 4.10 | 4.12^{58} |

S_{4} | 2A_{g} | (H → L)^{2} | 0.00 | 4.30 | 4.35 | 4.32 | … | … | 4.29^{58} |

S_{5} | 2B_{3u} | H → L + 1 | 1.73 | 4.18 | 4.35 | 4.43 | … | … | 4.66^{56} |

H-1 → L | |||||||||

S_{6} | 2B_{1g} | H-2 → L | 0.00 | 4.28 | 4.46 | 4.56 | 4.64 | 4.48 | 4.54^{58} |

H → L + 2 | |||||||||

S_{7} | 3B_{1g} | H-3 → L | 0.00 | 4.73 | 4.77 | 4.82 | 4.89 | 4.85 | 4.94^{58} |

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 2A_{g}, 1B_{1g}, and 1B_{3u} are already converged with respect to the active space size, the remaining states (in particular, both bright states 1B_{2u} and 2B_{3u}) 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 data^{56–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 D_{2h} symmetry where states of the same irreducible representation are far apart in energy and do not mix. Only in the case of the B_{1g} 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 S_{0} equilibrium geometry, all the excited states show a gradient only along the totally symmetric A_{g} 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, 8$\u2225$0, 0$\u2225$4, 8)/ANO-L-VDZP level^{39} [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 2B_{3u} and 2B_{1g} (Table II). Two-photon absorption experiments put the 2B_{1g} (4.54 eV) below the second bright state 2B_{3u} (4.66 eV) at the respective excited minimum. When the reorganization energies—predicted as ∼0.05 eV for 2B_{3u} and 0.23 eV for 2B_{1g} (Table III)—are considered, the state order is inverted in the Franck–Condon (FC) point.

. | S_{1}
. | S_{2}
. | S_{3}
. | S_{4}
. | S_{5}
. | S_{6}
. | ||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|

. | OPT^{39}
. | LVC . | OPT^{39}
. | LVC . | OPT . | LVC . | OPT . | LVC . | OPT^{39}
. | LVC . | OPT . | LVC . |

λ (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 |

### B. Wavefunction overlap calculations

Vibronic coupling between the considered diabatic states is observed both along totally symmetric A_{g} modes and along the symmetry-decreasing modes belonging to the B_{3u}, B_{2u}, and B_{1g} irreducible representations (Table I). As noted earlier, in D_{2h} 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 D_{2h} symmetry: this is particularly evident in the case of the first bright state S_{2}, which is the only B_{2u} state in D_{2h} 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 D_{2h} 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 D_{2h}-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.

. | . | Deviation from reference energy (eV) . | Absolute mean . | Standard . | ||||||
---|---|---|---|---|---|---|---|---|---|---|

Symmetry . | Level of theory . | S_{1}
. | S_{2}
. | S_{3}
. | S_{4}
. | S_{5}
. | S_{6}
. | S_{7}
. | deviation (eV) . | deviation (eV) . |

D_{2h} | SS(8:16) | 3.23 | 3.75 | 4.16 | 4.32 | 4.43 | 4.56 | 4.82 | … | … |

Irreducible rep. | A_{1} | B_{1} | B_{1} | A_{1} | A_{1} | B_{1} | B_{1} | |||

C_{2v} (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. | B_{2} | A_{1} | B_{2} | A_{1} | B_{2} | B_{2} | B_{2} | |||

C_{2v} (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. | B_{u} | B_{u} | A_{g} | A_{g} | B_{u} | A_{g} | A_{g} | |||

C_{2h} | 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 mean . | Standard . | ||||||
---|---|---|---|---|---|---|---|---|---|---|

Symmetry . | Level of theory . | S_{1}
. | S_{2}
. | S_{3}
. | S_{4}
. | S_{5}
. | S_{6}
. | S_{7}
. | deviation (eV) . | deviation (eV) . |

D_{2h} | SS(8:16) | 3.23 | 3.75 | 4.16 | 4.32 | 4.43 | 4.56 | 4.82 | … | … |

Irreducible rep. | A_{1} | B_{1} | B_{1} | A_{1} | A_{1} | B_{1} | B_{1} | |||

C_{2v} (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. | B_{2} | A_{1} | B_{2} | A_{1} | B_{2} | B_{2} | B_{2} | |||

C_{2v} (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. | B_{u} | B_{u} | A_{g} | A_{g} | B_{u} | A_{g} | A_{g} | |||

C_{2h} | 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 |

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 D_{2h} geometry: S_{4}/S_{5} along B_{3u} modes [$\Delta ESS(8:16)D2h$ = 0.11 eV], S_{5}/S_{6} along B_{2u} modes [$\Delta ESS(8:16)D2h$ = 0.13 eV], and S_{3}/S_{4} along B_{1g} modes [$\Delta 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 S_{6}–S_{5} (along B_{2u} modes) and S_{4}–S_{3} (along B_{1g} 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 S_{4} and S_{5} states along B_{3u} modes (i.e., A_{1} 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 D_{2h}-SS(8:16) value, and the S_{5}–S_{4} mixing is significantly reduced [0.012 at MS(8:16) and 0.006 at XMS(8:16) level, see Table V].

. | . | . | . | . | Deviation from . | $\u27e8Siref|Sjdispl\u27e9$ . | ||
---|---|---|---|---|---|---|---|---|

. | Modes . | Symmetry . | Level of theory . | ΔE (eV) . | reference ΔE (eV) . | Average . | Min . | Max . |

S_{5}–S_{4} | B_{3u} | C_{2v}(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 | |||

S_{6}–S_{5} | B_{2u} | C_{2v}(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 | |||

S_{4}–S_{3} | B_{1g} | C_{2h} | 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 from . | $\u27e8Siref|Sjdispl\u27e9$ . | ||
---|---|---|---|---|---|---|---|---|

. | Modes . | Symmetry . | Level of theory . | ΔE (eV) . | reference ΔE (eV) . | Average . | Min . | Max . |

S_{5}–S_{4} | B_{3u} | C_{2v}(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 | |||

S_{6}–S_{5} | B_{2u} | C_{2v}(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 | |||

S_{4}–S_{3} | B_{1g} | C_{2h} | 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 A_{1} representation along B_{3u} modes (C_{2v} symmetry), for which the bigger active space is needed to avoid artificially high S_{5}/S_{4} 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 B_{3u}-A_{1} states were computed with the bigger active space [i.e., MS(8:16) or XMS(8:16)].

### C. Accuracy of the LVC PES

The three different parameterizations of the LVC Hamiltonian will be named from now on LVC_{MS(16)}, LVC_{XMS(12)}, and LVC_{XMS(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 LVC_{MS(16)} diabatic PESs along A_{g} collective coordinates leading from the S_{0} minimum to the minima of the different LVC diabatic PESs (solid lines) with the energies of the corresponding adiabatic states recomputed at the D_{2h}-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 3B_{1g} and 2B_{3u} along the coordinate connecting the S_{0} and the 1B_{1g} 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.

To have a closer look at the performance of the LVC model in the minima, we consider the LVC_{MS(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 2A_{g} and 2B_{1g} and are 0.04 eV. At each minimum, also the energies of the other states are quite similar with the partial exceptions of states 2B_{3u} and 3B_{1g}, 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 D_{2h} symmetry. Note that since in D_{2h} 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 D_{2h} 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 2B_{3u}, 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 1B_{1g}/2B_{3u} crossing actually corresponds to a quasi-triple CoI involving also the 2A_{g} state and reproduces the absolute energies up to 0.02 eV. A further quasi-triple CoI involving the 1B_{3u}, 1B_{2u}, and 1B_{1g} states (proposed previously based on orbital analysis and CoI search^{39}) 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 1B_{3u}/2A_{g} crossing, four states are found in <0.17 eV (i.e., also 1B_{2u} and 1B_{1g}), suggesting that a quasi-fourfold CoI might exist in the proximity of that structure.

. | . | States . | ||||||
---|---|---|---|---|---|---|---|---|

. | . | S_{1}
. | S_{2}
. | S_{3}
. | S_{4}
. | S_{5}
. | S_{6}
. | S_{7}
. |

CoI . | Methods . | 1B_{3u}
. | 1B_{2u}
. | 1B_{1g}
. | 2A_{g}
. | 2B_{3u}
. | 2B_{1g}
. | 3B_{1g}
. |

1B_{3u}/1B_{2u} | 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 | |

1B_{3u}/1B_{1g} | 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 | |

1B_{2u}/1B_{1g} | 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 | |

1B_{1g}/2B_{3u} | 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 | |

2A_{g}/2B_{3u} | 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 | |

2B_{3u}/2B_{1g} | 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 . | ||||||
---|---|---|---|---|---|---|---|---|

. | . | S_{1}
. | S_{2}
. | S_{3}
. | S_{4}
. | S_{5}
. | S_{6}
. | S_{7}
. |

CoI . | Methods . | 1B_{3u}
. | 1B_{2u}
. | 1B_{1g}
. | 2A_{g}
. | 2B_{3u}
. | 2B_{1g}
. | 3B_{1g}
. |

1B_{3u}/1B_{2u} | 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 | |

1B_{3u}/1B_{1g} | 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 | |

1B_{2u}/1B_{1g} | 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 | |

1B_{1g}/2B_{3u} | 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 | |

2A_{g}/2B_{3u} | 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 | |

2B_{3u}/2B_{1g} | 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 |

### D. Dynamics of electronic populations

Figure 4 shows the time evolution of the electronic populations up to 2 ps after the initial photo-excitation to either the first (1B_{2u}) or the second (2B_{3u}) bright states according to the LVC_{MS(16)} and LVC_{XMS(16)} parameterizations (results with LVC_{XMS(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. LVC_{MS(16)} and LVC_{XMS(16)} Hamiltonians deliver similar predictions: 1B_{2u} decays essentially on the lowest state 1B_{3u}, while after an initial excitation to 2B_{3u}, 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 1B_{2u}, which reaches its maximum population (∼0.5) in 100 fs and then slowly decays toward 1B_{3u}. The intermediate population of 1B_{2u} is consistent with the two-step interpretation of Borrego-Varillas *et al.* who reported transient signatures of 1B_{2u} 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 2B_{3u} agrees with experimental time constants reported in the literature.^{36,39,40}

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

. | MS(8:16) . | XMS(8:16) . | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

State . | 1B_{3u}
. | 1B_{2u}
. | 1B_{1g}
. | 2A_{g}
. | 2B_{3u}
. | 2B_{1g}
. | 3B_{1g}
. | 1B_{3u}
. | 1B_{2u}
. | 1B_{1g}
. | 2A_{g}
. | 2B_{3u}
. | 2B_{1g}
. | 3B_{1g}
. |

1B_{3u} | 0.159 | 0.159 | ||||||||||||

1B_{2u} | 0.043 | 0.184 | 0.030 | 0.184 | ||||||||||

1B_{1g} | 0.199 | 0.196 | 0.257 | 0.111 | 0.096 | 0.257 | ||||||||

2A_{g} | 0.108 | 0.116 | 0.049 | 0.257 | 0.105 | 0.043 | 0.058 | 0.256 | ||||||

2B_{3u} | 0.108 | 0.124 | 0.027 | 0.054 | 0.126 | 0.072 | 0.176 | 0.096 | 0.028 | 0.126 | ||||

2B_{1g} | 0.087 | 0.096 | 0.037 | 0.235 | 0.152 | 0.266 | 0.056 | 0.126 | 0.059 | 0.109 | 0.146 | 0.267 | ||

3B_{1g} | 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) . | ||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

State . | 1B_{3u}
. | 1B_{2u}
. | 1B_{1g}
. | 2A_{g}
. | 2B_{3u}
. | 2B_{1g}
. | 3B_{1g}
. | 1B_{3u}
. | 1B_{2u}
. | 1B_{1g}
. | 2A_{g}
. | 2B_{3u}
. | 2B_{1g}
. | 3B_{1g}
. |

1B_{3u} | 0.159 | 0.159 | ||||||||||||

1B_{2u} | 0.043 | 0.184 | 0.030 | 0.184 | ||||||||||

1B_{1g} | 0.199 | 0.196 | 0.257 | 0.111 | 0.096 | 0.257 | ||||||||

2A_{g} | 0.108 | 0.116 | 0.049 | 0.257 | 0.105 | 0.043 | 0.058 | 0.256 | ||||||

2B_{3u} | 0.108 | 0.124 | 0.027 | 0.054 | 0.126 | 0.072 | 0.176 | 0.096 | 0.028 | 0.126 | ||||

2B_{1g} | 0.087 | 0.096 | 0.037 | 0.235 | 0.152 | 0.266 | 0.056 | 0.126 | 0.059 | 0.109 | 0.146 | 0.267 | ||

3B_{1g} | 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 2B_{3u}, the initial decay (∼10 fs) is toward 2B_{1g} and 2A_{g} according to LVC_{MS(16)} and toward 2B_{1g}, 1B_{1g}, and directly 1B_{2u} according to LVC_{XMS(16)}. These differences can be attributed to corresponding differences in the pattern of the couplings reported in Table VII. Indeed, the couplings of 2B_{3u} with 1B_{1g} and 1B_{2u} are remarkably larger according to LVC_{XMS(16)}. On the contrary, the coupling of 2B_{3u} with 2A_{g} is larger according to LVC_{MS(16)}. The latter also predicts a much larger coupling of the higher-energy state 2B_{1g} with 2A_{g} explaining why, despite its energy, 2B_{1g} gains some transient population, which, according to LVC_{MS(16)}, reaches slightly larger values and decays at a slightly slower rate than in the case of LVC_{XMS(16)}.

Analysis of Fig. 4 suggests that after photoexcitation to 1B_{2u} the dynamics is quite simple, being essentially characterized by a progressive (approximatively mono-exponential) flow of population from 1B_{2u} to the lowest-energy state 1B_{3u}. This is not surprising considering that at the FC position, the third state, 1B_{1g} is ∼0.5 eV higher in energy than 1B_{2u}. However, Table VII shows that 1B_{1g} is strongly coupled to both 1B_{2u} and 1B_{3u} states. More specifically, the norm of its coupling to these two states is, respectively, more than three [LVC_{XMS(16)}] and more than four [LVC_{MS(16)}] times larger than the direct 1B_{1g}/1B_{2u} coupling. A small transient population on 1B_{1g} is actually seen in Fig. 4 for the Hamiltonian with the larger couplings [LVC_{MS(16)}]. In Fig. 5, we investigate in greater detail the impact on the 1B_{2u} → 1B_{3u} transfer of the existence of 1B_{1g} 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 “1B_{2u} + 1B_{3u},” the three-state model “1B_{2u} + 1B_{3u} + 1B_{1g},” and the six-state model obtained including all states except 1B_{1g}. Differences are striking: according to the two-state model, the population transfer is much slower, smaller in amplitude, and shows large oscillations. Including also 1B_{1g}, 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 1B_{1g} is removed. In the six-state-model, the predicted population flow from 1B_{2u} to 1B_{3u} is, in fact, similar to what is obtained with the complete seven-state model. Actually, in the long-time limit, 1B_{3u} 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 1B_{1g} has a dramatic impact on the 1B_{2u} → 1B_{3u} 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 1B_{1g} is partially contrasted by the higher-energy states, which slow down the rise of the population of 1B_{3u}. On the long-time scale, however, according to the seven-state model, 1B_{1g} maintains a weak population (∼3%). If such a state is not included in the calculation, this small population flows to 1B_{3u}, making the yield of this state even larger (six-state model). 1B_{1g} and higher-energy states play a qualitatively similar role also according to the XMS(8:16) parameterization, but couplings with 1B_{1g} are smaller. In conclusion, the faster 1B_{2u} → 1B_{3u} decay predicted by LVC_{MS(16)} with respect to LVC_{XMS(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 1B_{1g} (see Table VII).

Figure 6 plots the diabatic LVC PES at the average position of the WP as a function of time according to the LVC_{MS(16)} Hamiltonian [results for LVC_{XMS(16)} are very similar and are given in Fig. S12 of the supplementary material]. It shows that at all times, S_{1} and S_{2} are well separated in energy and rather distant from two pairs of close-lying states, namely, S_{3}–S_{4}, and S_{5}–S_{6}. 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 1B_{2u}, since the initial potential energy of the WP is 3.75 eV (Table II) and the lowest 1B_{1g}/1B_{2u} crossing is at ∼4.2 eV (Table VI), it is noteworthy that the same picture applies also for an initial excitation to 2B_{3u} although several crossings between diabatic states are reachable at this energy, including the (quasi) triple-crossings 1B_{3u}/1B_{2u}/1B_{1g} and 1B_{1g}/2A_{g}/2B_{3u}.

Finally, Fig. 7 reports the expectation values of all the total-symmetric modes as a function of time for the LVC_{MS(16)} Hamiltonian. The results for LVC_{XMS(16)} are shown in Fig. S13 and are very similar. Both starting from 1B_{2u} and 2B_{3u}, 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 1B_{2u} and 2B_{3u},^{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 A_{g} 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.

We conclude this section mentioning that LVC_{XMS(12)} predicts a very different dynamics (Fig. S14 in the supplementary material), characterized by the fact that both starting from 1B_{2u} and 2B_{3u}, the states 2B_{3u} and 2A_{g} 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).

## V. CONCLUSIONS

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 (2A_{g}) 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 2A_{g} 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 approach^{73} 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., 2A_{g}, already well described with the full-*π* active space, and the 2B_{3u} 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. LVC_{XMS(16)} and LVC_{MS(16)} dynamics are qualitatively similar. Still, both for an initial excitation to 1B_{2u} and to 2B_{3u}, LVC_{MS(16)} predicts that the decay from 1B_{2u} to 1B_{3u} 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 1B_{2u} (S_{2}), the population progressively flows to 1B_{3u} (S_{1}). In particular, the population growth with a sub-100 fs time constant predicted by the LVC_{MS(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 1B_{1g}, 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 2B_{3u} (S_{5}) leads to its ultrafast (sub-100 fs) depopulation in favor of a number of intermediate states, especially 1B_{2u}, followed by a much slower progressively decay to 1B_{3u} (S_{1}), 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 1B_{2u} (S_{2}) → 1B_{3u} (S_{1}) 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}

## SUPPLEMENTARY MATERIAL

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 LVC_{MS(16)} and LVC_{XMS(16)} Hamiltonians, average position of the A_{g} modes during wavepacket propagations for the LVC_{XMS(16)} Hamiltonian, comparison of the population dynamics with the LVC_{MS(16)}, LVC_{XMS(12)}, and LVC_{XMS(16)} Hamiltonians, and norm of the diabatic coupling vectors for LVC_{XMS(12)} parameterization.

## AUTHORS’ CONTRIBUTIONS

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

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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

## REFERENCES

^{*}/nπ

^{*}internal conversion in thymine

^{*}→ nπ

^{*}decay in 9H-adenine. A quadratic vibronic coupling model

^{*}/nπ

^{*}decay of thymine in water as a test case

^{*}/nπ

^{*}decay of the epigenetic nucleobase 1,5-dimethyl-cytosine in the gas phase

^{*}/nπ

^{*}population transfer in uracil and 5-fluoro-uracil in water and acetonitrile

^{1}L

_{a}to

^{1}L

_{b}in pyrene: A theoretical study

_{2}

^{*}states in the photophysics of pyrazine: A quantum dynamics study