Electronic Structure Methods for Simulating Flavin’s Spectroscopy and Photophysics: Comparison of Multi-reference, TD-DFT, and Single-Reference Wave Function Methods

The use of flavins and flavoproteins in photocatalytic, sensing, and biotechnological applications has led to a growing interest in computationally modeling the excited-state electronic structure and photophysics of flavin. However, there is limited consensus regarding which computational methods are appropriate for modeling flavin’s photophysics. We compare the energies of low-lying excited states of flavin computed with time-dependent density functional theory (TD-DFT), equation-of-motion coupled cluster (EOM-EE-CCSD), scaled opposite-spin configuration interaction [SOS-CIS(D)], multiconfiguration pair-density functional theory (MC-PDFT), and several multireference perturbation theory (MR-PT2) methods. In the first part, we focus on excitation energies of the first singlet excited state (S1) of five different redox and protonation states of flavin, with the goal of finding a suitable active space for MR-PT2 calculations. In the second part, we construct two sets of one-dimensional potential energy surfaces connecting the S0 and S1 equilibrium geometries (S0–S1 path) and the S1 (π,π*) and S2 (n,π*) equilibrium geometries (S1–S2 path). The first path therefore follows a Franck–Condon active mode of flavin while the second path maps crossings points between low-lying singlet and triplet states in flavin. We discuss the similarities and differences in the TD-DFT, EOM-EE-CCSD, SOS-CIS(D), MC-PDFT and MR-PT2 energy profiles along these paths. We find that (TD-)DFT methods are suitable for applications such as simulating the spectra of flavins but are inconsistent with several other methods when used for some geometry optimizations and when describing the energetics of dark (n,π*) states. MR-PT2 methods show promise for the simulation of flavin’s low-lying excited states, but the selection of orbitals for the active space and the number of roots used for state averaging must be done carefully to avoid artifacts. Some properties, such as the intersystem crossing geometry and energy between the S1 (π,π*) and T2 (n,π*) states, may require additional benchmarking before they can be determined quantitatively.


■ INTRODUCTION
Flavins are cofactors in several classes of enzymes and in photoreceptor proteins, where they are typically bound as either flavin mononucleotide (FMN) or flavin adenine dinucleotide (FAD).Flavins are often encountered in one of three redox states: the oxidized quinone (Fl), the one-electron reduced radical semiquinone (FlH • ), and the two-electron fully reduced hydroquinone (FlH 2 ).−3 Therefore, the semiquinone and hydroquinone redox states are frequently also encountered in their deprotonated anionic forms, Fl •− and FlH − , under physiological conditions (Figure 1).These redox and protonation states play an important role in either coupled or stepwise electron and/or proton transfer reactions (Figure 1).Here, we refer collectively to FMN, FAD, riboflavin, and lumiflavin (LF) as flavins, since they all share the same tricyclic 7,8-dimethyl alloxazine group.
−7 Flavoproteins are critical to all the main energy generation metabolisms 8 including photosynthesis, 9 aerobic and anaerobic respiration, 10−12 and denitrification. 13,14Flavins also serve as chromophores in several photoreceptors and light-responsive enzymes including cryptochromes, 15 LOV domains (phototropins), 16 BLUF, 17 and DNA photolyases. 18,19n its oxidized form, flavin has two absorption bands in the near-UV and visible range, one band centered around 450 nm, associated with excitation to the first singlet excited state (S 1 ), and a second band at around 380 nm associated with excitation to a higher singlet excited state.These absorption wavelengths are relatively insensitive to the solvent or protein environment. 20,21Both of those states have π,π* excitation character.Other low-lying excited states also exist in that energy range, including triplet states and singlet states, but those states have substantially smaller oscillator strengths (extinction coefficients) and do not appear in the UV/visible spectrum.Even though such dark states do not get directly excited, they can play an important role in flavins' photophysics, as exemplified by LOV domains that undergo intersystem crossing (ISC) followed by adduct formation. 22,23lue light absorption excites flavin to its S 1 state.Higher energy (near-UV or UV) absorption also ultimately populates the same state after a series of internal conversions and vibrational relaxations, per Kasha's rule. 24Therefore, flavin's S 1 state often acts as the launch pad for flavin's further photophysics.−35 Unlike the π,π* transitions responsible for flavins' two UV/ visible absorption bands, the energies of flavin's low-lying dark states, especially those having n,π* character, can be sensitive to interactions with a solvent or protein. 21,36,37Since those states are relevant to flavin's photophysics, computational modeling can be a powerful tool to study their energetics in different protein environments.However, most computational studies on flavins have focused on their ground state or on bright excited states.−39 On the other hand, TD-DFT has not been as thoroughly tested for describing excited states relevant to flavin's photophysics, including triplet, n,π*, or redox states.
When dealing with a manifold of excited states that are close together in energy, multiconfigurational wave function methods such as complete active-space self-consistent field (CASSCF) 40 are suitable.In CASSCF, a full configuration interaction is carried out for select electrons and orbitals (the active space).While CASSCF misses dynamical electron correlation for electrons outside of the active space, post-CASSCF (multireference) methods can be used to correct for this missing correlation.Second-order multireference perturbation (MR-PT2) theories such as single state (SS-)CASPT2 and multistate (MS-)CASPT2 are widely used. 41,42More recently, to address issues with artifacts in SS-CASPT2 and MS-CASPT2 arising at near-degeneracies of the underlying CASSCF reference wave function, methods like the extended multistate (XMS) CASPT2 43 and extended dynamically weighted (XDW) CASPT2 44 have been developed.The former employs a treatment of the zeroth order Hamiltonian proposed as by Granovsky 45 which removes the artifacts but alters the energy of the states of interest compared to MS-CASPT2.The latter approach, XDW-CASPT2, also addresses the artifacts but largely preserves the energies of the MS-CASPT2 states.
An alternative approach for accounting for static and dynamical electron correlation is the multiconfiguration pairdensity functional theory (MC-PDFT), 46,47 which combines favorable characteristics of CASSCF and DFT.MC-PDFT has been used in a range of applications from transition metal chemistry to mechanistic studies on biological photoreceptors. 48,49−39,50−81 For multireference calculations, the issue of active space selection is a recurring theme and a decision is often made that balances computational cost and computa- The Journal of Physical Chemistry B tional accuracy due to the steep increase in computing time with the size of the active space.−86  and Hammes-Schiffer applied their automated protocol that is based on Huckel theory for the selection of a suitable active space for flavin's π-conjugated orbitals. 74However, the choice of active space is also dependent on the problems and states of interest.For example, for calculating bright π,π* transitions, it is often adequate to include only π electrons and π and π* orbitals in the active space.For the calculation of dark n,π* states, however, including nonbonding orbitals in the active space is required.Therefore, there is still a need to tailor active spaces not only to the system but also to the problem of interest.
We present a benchmark study focused on flavin's low-lying excited states that are relevant to flavin's early photophysics.We focus on LF (Figure 1), a minimal flavin model that contains the spectroscopically and photochemically relevant isoalloxazine tricyclic moiety.We test several TD-DFT functionals, equation-of-motion coupled cluster (EOM-EE-CCSD), 87,88 scaled opposite-spin configuration interaction [SOS-CIS(D)], 89 MC-PDFT, and several MR-PT2 methods.For MR-PT2, we present the effect of modifying parameters including the active spaces and state averaging.We initially focus on excitation energies to the S 1 excited state of each redox state.We then move on to potential energy scans (PESs) of low-lying excited states, tracing the coordinates from the Franck−Condon point (at the ground state optimized geometry) to the S 1 (π,π*) minimum and then from the S 1 (π,π*) minimum to the S 2 (n,π*) minimum.These scans are used to probe the reliability of different methods across key geometries and state crossing points.
■ METHODS Vertical Excitation Calculations.Each of the five redox states of LF shown in Figure 1 (Fl, FlH • , Fl •− , FlH − , FlH 2 ) were optimized in the gas phase using second-order Møller− Plesset perturbation theory (MP2) with the correlationconsistent polarized valence triple-ζ (cc-pVTZ) basis set.The optimizations were performed without the use of symmetry.In the case of Fl, FlH • , and Fl •− , the isoalloxazine ring remains planar during the optimization.−92 Vertical excitation energies (VEEs) were computed using CASSCF, MC-PDFT, SS-CASPT2, MS-CASPT2, XMS-CASPT2, and XDW-CASPT2.Energies were also computed using TD-DFT with the B3LYP functional, which has been shown to give adiabatic excitation energies (AEEs) for the oxidized form of flavin consistent with experiments. 21,38,39,61he TD-DFT calculations were carried out using cc-pVTZ, while the ANO-L-VDZP basis set was used for the CASSCF, MC-PDFT, and MR-PT2 calculations.
To determine a suitable starting point for the active space, orbital entanglement calculations were carried out at the CASSCF density matrix renormalization group (CASSCF-DMRG) level of theory 86,93 and ANO-L-VDZP basis set.Orbital entanglement figures were generated using autoCAS, 94 and served as starting point for the selection of the active space orbitals for CASSCF and multireference calculations.The number of active space orbitals and electrons were gradually reduced based on their occupancy to test the effect of using a smaller active space on the excitation energies.The number of roots used in the state averaging was also benchmarked.
To mitigate the issue of intruder states in MR-PT2 calculations, an imaginary shift of 0.2 hartree was used. 95We tried using a real shift, which gave similar energies except for a The Journal of Physical Chemistry B few spurious spikes in the excitation energy for some geometries. 96,97We also tested the effect of applying a shifted zeroth-order Hamiltonian (IPEA shift) utilizing the default value of 0.25 atomic units. 98The IPEA shift is often used as an approximate correction for an unbalanced description of openshell and closed shell electronic states, although it may not be needed if a sufficiently large basis set and active space is used. 99otential Energy Scans.To generate one-dimensional PESs along coordinates relevant to flavin's photophysics, we reoptimized the S 0 , S 1 , and S 2 states of flavin using (TD-)B3LYP/cc-pVTZ with C s symmetry enforced.We then generated two 1-D PESs by linear interpolation of Cartesian coordinates.The first PES (S 0 −S 1 path, see Figure 2 left) connects the S 0 optimized geometry to the S 1 optimized geometry and therefore follows a Franck−Condon-active mode.The second PES (S 1 −S 2 path, see Figure 2 right) connects the S 1 optimized geometry to the S 2 optimized geometry and therefore traces the path where the two states may potentially cross.We also added additional points on either side of each path through a linear extrapolation.In total, each path contains twenty-one geometries obtained through interpolation and extrapolation, as shown in Figure 2.
Single-point calculations along the PESs were performed using TD-DFT with the B3LYP, 100,101 camB3LYP, 102 and ωB97X-D 103 functionals and the cc-pVTZ basis set.EOM-EE-CCSD and SOS-CIS(D) energies were computed with the cc-pVDZ basis set.MR-PT2 calculations along the scans were also carried out using the ANO-L-VDZP basis set.A [14,12]  active space was used for the PESs, comprising two nonbonding orbitals with electron density on the N 1 and N 5 nitrogen atoms and five π and five π* orbitals.3-Root state averaging was used for each symmetry and spin.
MP2 and DFT calculations were run using Gaussian 16. 104 Multiconfigurational calculations were performed with Open-Molcas version 22.10. 105The density matrix renormalization group (DMRG) calculations for orbital entanglement were run using the QCMaquis interface with OpenMolcas. 106The resultant orbital entanglement data were analyzed using autoCAS. 94EOM-EE-CCSD and SOS-CIS(D) calculations were run using Q-Chem 5.4. 107

■ RESULTS AND DISCUSSION
The geometries for Fl, FlH • , Fl •− , FlH − , FlH 2 optimized at the MP2/cc-pVTZ level of theory are compared to the geometries optimized at the B3LYP/cc-pVTZ level of theory from ref 21.In Figure 3, we align the molecule along the x−y plane as closely as possible and find the maximum and minimum displacements of heavy atoms along the z-axis.This distance, labeled d, indicates the extent of the "butterfly" bending of the isoalloxazine rings in the MP2 and B3LYP geometries for the different redox states.The MP2 and B3LYP geometries are also shown superimposed in the same figure.Even though all molecules were optimized without symmetry constraints, Fl, FlH • , and Fl •− remain planar, with only minimal deviations from planarity.The MP2 and B3LYP geometries are similar, which is also reflected in the TD-B3LYP S 1 excitation energies that differ by no more than 0.05 eV for the two geometries.However, FlH − and FlH 2 show significantly more bending with MP2 compared to B3LYP.The distance d is 1.60 Å for FlH − at the MP2 level of theory compared to 1.25 Å with B3LYP, while in FlH 2 , d is 1.32 Å with MP2 compared to 0.76 Å with B3LYP.This causes a larger difference in the excitation energies computed with the MP2 and B3LYP geometries (a 0.22 eV difference).Therefore, accurate calculations of thermodynamic quantities or excitation energies for those two reduced states of flavin may necessitate the use of geometries optimized with MP2 or a better level of theory.
Figure 4 presents the orbital entanglement diagram for the oxidized form of flavin, which was generated from CASSCF-DMRG wave functions with a [38,36] active space.The diagrams feature a circular arrangement of active orbitals where strong orbital "entropy," a term used in this context to signify a high degree of correlation with other orbitals, is indicated by larger gray circles and mutual information is indicated by the line thickness between the orbitals.A set of 15 orbitals, highlighted in the red boxes in the diagram, exhibited higher single-orbital entropies and mutual information.These orbitals comprise eight bonding π orbitals (16 electrons) and seven antibonding π* orbitals, which are shown in Figure S1a of the Supporting Information document.Consequently, we employed CASSCF with a 16 electron and 15 orbital active space as the starting point for further benchmarking.
While the DMRG-CASSCF calculations indicate the relevance of these 15 orbitals to the molecule's multiconfigura- The Journal of Physical Chemistry B tional electronic structure, such a large active space makes routine calculations on flavins intractable.This is especially true if one is interested in looking at other excited states such as n,π*, which would necessitate the addition of n orbitals on top of π and π*.Therefore, we examined the effect of reducing the active space on the lowest excitation energy of flavin.
In Figure 5, MR-PT2 and MC-PDFT excitation energies are plotted as a function of active space.3-Root state averaging is used for the underlying CASSCF wave function in all cases.The effect of using additional roots is explored further in Supporting Information Figure S2, where it is shown that the MR-PT2 energies are relatively stable with the use of three or more states in the CASSCF state averaging.We started with a 16 electron and 15 orbital [16,15] active space shown in Figure S1a.Removing orbitals based on their occupancy had a relatively small impact on S 0 −S 1 excitation energies down to a [10,10] active space.However, removing additional orbitals resulted in significant energy variations, indicating unbalanced electron correlation in smaller spaces.A minimal active space of [2,2] gives reasonable excitation energies for most MR-PT2 methods, as has been also shown by Udvarhelyi and Domratcheva. 53,57Overall, the results in Figure 5 suggest that a [10,10] active space is sufficient for accurate π,π* excitation energy calculations when using the orbitals labeled as π 1 , π 2 , π 3 , π 4 , π 5 , π 1 *, π 2 *, π 3 *, π 4 *, and π 5 * in Figure S1a.SS-CASPT2, MS-CASPT2, and XDW-CASPT2 provide largely consistent energies.XMS-CASPT2 gives energies that are consistently 0.2 eV above those methods, while MC-PDFT gives energies that are around 0.1 to 0.25 eV below these methods, depending on the active space used.
In Figure 6, we show AutoCAS maps constructed from CASSCF-DMRG wave function calculations on FlH • , Fl •− , FlH 2 , and FlH − .The calculations for FlH • and Fl •− use a [39,36] active space, while those for FlH 2 and FlH − use a [40,36] active space.We find that eight bonding and six antibonding orbitals are important for both FlH • and Fl •− .For FlH 2 and FlH − , nine bonding and six antibonding orbitals with high entropy were found instead.These orbitals are highlighted in red boxes and are shown in Supporting Information Figure S1b−e along with their occupancies.
We computed the excitation energies for these four redox and protonation states as a function of decreasing active space size by sequentially removing the highest occupancy orbitals, as we did for Fl.The excitation energies obtained from these calculations, along with TD-B3LYP energies, are shown in Figure 7.
For FlH • and Fl •− , the MR-PT2 and MC-PDFT calculations with large active spaces gave results largely similar to TD-    [39,36] actives spaces were used for the semiquinones and [40,36] for the hydroquinones.Gray circle size represents the degree of orbital entropy while line thickness indicates mutual information between orbital pairs.Red boxes highlight orbitals with high entropy.
The Journal of Physical Chemistry B B3LYP.For FlH • , reducing the active space size had a limited effect down to a [7,7] active space.A minimal active space, [3,3], also gives similar results, but a [5,5] active space appears imbalanced.In FlH − , on the other hand, any active space [3,3]  or larger gave similar energies.
In FlH 2 , and FlH − , which have nonplanar geometries, the convergence of the excitation energies with increasing active space size is much less smooth compared to the other states.Up to the largest active space tested, [16,14], there continue to be changes in the excitation energies.This may be attributed to the bent structure of these two redox states, which allow mixing between nonbonded or σ orbitals with π orbitals that likely results in stronger correlation between a larger set of orbitals.This can be seen to some extent in the autoCAS maps in Figure 6, where FlH 2 and FlH − have more orbitals involved in the mutual information exchange, even though the degree of entanglement is smaller for each orbital.
In FlH 2 , there is a large difference of ca.1.0 eV between the MC-PDFT and SS-CASPT2 results compared to those calculated with MS/XMS/XDW-CASPT2.Such differences between SS-CASPT2 and MS-CASPT2 have been previously discussed 45 and are usually resolved through the XMS or XDW extension to these theories, as appears to be the case here.
PESs for Flavin's Low-Lying Excited States.For the application of electronic structure methods toward studying the photophysics of flavin, it is important to move beyond single point energy calculations and to explore the energies of different methods along relevant coordinates.Here, we focus on PESs near flavin's Franck−Condon region along two modes (Figure 2).One mode (S 0 −S 1 ) connects the geometries of the S 0 and S 1 minima optimized at the (TD)-B3LYP/cc-pVTZ level of theory.The second mode (S 1 −S 2 ) connects the S 1 and S 2 minima optimized at the TD-B3LYP/cc-pVTZ level of theory.The aim is to check the sensitivity of the relative energies, curvatures, and crossing points of the low-lying singlet and triplet excited states along these paths to the electronic structure method.
In Figure 8, we show MS-CASPT2 energies for five states (S 0 , S 1 , S 2 , T 1 , and T 2 ) along the two paths.The [14,12] active space comprises the 5 π and 5 π* orbitals consistent with the [10,10] active space from Figure 5 plus two nonbonding orbitals localized on the N 1 and N 5 nitrogen atoms.S 1 is the optically active π,π* state, while S 2 is the optically inactive n,π* state involving the nonbonding orbitals.The triplet states, T 1 and T 2 , have the same electronic transitions as the S 1 and S 2 states, respectively.The S 1 state crosses with both the S 2 and the T 2 states along the S 1 −S 2 path.
We first focus on the S 0 and S 1 PESs along the S 0 −S 1 path, since those two states are responsible for the spectroscopic properties of flavin's first excitation band.We then shift our focus to the S 1 , S 2 , and T 2 PESs along the S 1 −S 2 path, since their crossings represent possible internal conversion or ISC points, respectively.
Figure 9 plots the S 0 and S 1 energies along the S 0 −S 1 path with different methods.For S 0 , CASSCF stands out from other  The Journal of Physical Chemistry B curves at large distortion.This is likely due to missing dynamical electron correlation, particularly σ orbitals that become involved with bond stretching.However, even among other methods that do account for dynamical electron correlation, we find that there are differences in the positions of the S 0 minima along the PES.CCSD is the only other method that has the same minimum point as the B3LYPoptimized S 0 state.Conversely, camB3LYP and ωB97X-D give geometries that are further away from the S 1 minimum and fall along the extrapolated portion of the PES.MR-PT2, MC-PDFT, and SOS-CIS(D) methods, on the other hand give minima that are closer to the S 1 equilibrium structure.Nonparallelity errors (NPEs), which are reported relative to the XDW-CASPT2 (IPEA = 0.25) energy profile, reflect those differences as well (Table 1).NPEs are calculated as the largest deviation from the reference XDW-CASPT2 curve minus the smallest deviation from the same curve.Methods that have similar minimum points to XDW-CASPT2 and have similar curvature have small S 0 NPEs, while methods with different minimum energy geometries display larger S 0 NPEs.
The differences between methods become substantially more pronounced when looking at the S 1 state, both in terms of the excitation energy as well as the shape (curvature) of the excited state potential energy surface along this path.VEEs, which are indicated with circles in the top left section of Figure 9, reflect the S 0 −S 1 energy differences at the same position where a method's S 0 minimum is found.Even if we exclude CASSCF from the analysis due to its missing dynamical electron correlation, we find that VEEs vary from 2.59 eV (MS-CASPT2 with IPEA = 0) to 3.49 eV (EOM-EE-CCSD), a range of 0.90 eV.Changes in the curvature of the excited state potential energy surface are reflected by the S 1 NPEs shown in Table 1.While some excited state methods have PESs that resemble the XDW-CASPT2 PES (S 1 NPE within 0.11 eV), several methods stand out has having a different curvature along this path.Specifically, SS-CASPT2 (NPE = 0.37−0.42eV), SOS-CIS(D) (NPE = 0.41 eV), MC-PDFT (NPE = 0.28 eV), and XMS-CASPT2 (NPE = 0.19 eV) all have a larger curvature compared to XDW-CASPT2.On the other hand.TD-DFT methods, which give NPEs ranging from 0.10 to 0.18, have smaller curvature.
In the top right part of Figure 9, circles are used to indicate the S 1 minima for the different methods calculated along the S 0 −S 1 path.The ΔE value at those points reflects the AEEs for the different methods.We note that those are approximate AEEs because they are not true minima but instead are the lowest energy points computed along a TD-B3LYP PES.Nonetheless, there is a marked variation in these energies, which is a result of both the differences in the VEEs as well as the differences in curvatures of the different methods.AEEs along this path range from 2.33 eV (CASPT2 with IPEA = 0) to 3.19 eV for EOM-EE-CCSD, a 0.86 eV range.
In Figure 10, we plot S 0 , S 1 , S 2 , and T 2 energies along the S 1 −S 2 path.The circles on the S 0 PES again indicate differences in the positions of the minima for the different methods, which all lie somewhere between the S 1 and S 2 geometries.In the top left panel of Figure 10, we show both the S 1 and S 2 states.For all methods shown, the two states cross at some point along the S 1 −S 2 path.In all cases, the S 1 minimum remains lower in energy than the S 2 minimum, suggesting (at least, in the gas phase) that the S 2 (n,π*) state would not get populated from the S 1 state.Since solvents typically stabilize (red-shift) π,π* excited states and destabilize (blue-shift) n,π* states, 108 we expect the same conclusion will also apply to flavins in solution.The S 1 −S 2 crossings for all methods, indicated using circles, are geometrically closer to the S 2 minima and often lie well above the S 1 minima.In Table 2 we tabulate ΔE Sd 1 /Sd 2 , which are energy differences between the S 1 −S 2 crossing points and the S 1 minima, for the different methods.
The S 1 /T 2 plots (Figure 10 right) tell a different story.There, the relative energies of the S 1 and T 2 minima are more comparable, with some methods showing a lower S 1 state minimum energy point and other methods showing a lower T 2 minimum.For example, TD-DFT methods show that the T 2 minimum is lower than the S 1 minimum, while MR-PT2 methods show that the S 1 minimum is lower than T 2 .The S 1 /  The Journal of Physical Chemistry B T 2 crossings lie somewhere in the middle between the two coordinates.As discussed by Salzmann et al. in multiple studies, 36,76,77 solvent effects or a protein environment can play an important role in changing flavin's photophysics by modulating the relative energies of the S 1 and T 2 states.However, here we also find that the relative energies of those states are sensitive to the choice of computational method.Table 2 presents ΔE Sd 1 /Td 2 , the relative energies of the S 1 /T 2 ISC point compared to the S 1 minimum along the S 1 −S 2 path, alongside the corresponding ΔE Sd 1 /Sd 2 values.We note that these crossing points were obtained from the interpolation and do not necessarily represent the minimum energy crossing points.However, the variations in these numbers for the different methods indicates either that (1) the barrier to those crossings is different for the different methods, or (2) the minimum energy crossing geometry varies from one method to another.For example, along the S 1 −S 2 path, TD-DFT methods indicate that the S 1 /T 2 ISC point coincides with the S 1 minimum geometry (with a barrier of ∼0 eV) while SOS-CIS(D) gives a ΔE Sd 1 /Td 2 of 0.62 eV.EOM-EE-CCSD gives a ΔE Sd 1 /Td of 0.10 eV while MR-PT2 methods give ΔE Sd 1 /Td 2 values that vary from 0.23 to 0.36 eV.
When looking more generally at the trends in both ΔE Sd 1 /Sd 2 and ΔE Sd 1 /Td 2 , most MR-PT2 methods and MC-PDFT appear to treat the π,π* and n,π* states on a similar footing, since the crossing geometries do not change considerably from one method to the other.The IPEA shift has a large effect on the S 1 , S 2 , and T 2 energies relative to the reference S 0 state, but it does not have much of an impact on their energies relative to each other since all three states have similarly open-shell character.This is results in ΔE Sd 1 /Sd 2 and ΔE Sd 1 /Td 2 barriers for the IPEA = 0 calculations that are comparable in magnitude to their IPEA = 0.25 counterparts.TD-DFT methods, on the other hand, predict lower energies for both the S In Figure 11, we show the energy profiles computed using MS-CASPT2 (IPEA = 0.25) with different state averaging or active space along the S 0 −S 1 path.Using 3 or more states in   A more comprehensive list of computational studies can be found in a recent review article by Kar, Miller, and Mroginski. 78These studies employ a range of computational methods including semiempirical, TD-DFT, coupled cluster, and multiconfigurational approaches.Even among multiconfigurational studies, methodological details such as the number of roots used in state averaging or the size of the active space may differ.With the continued advancement in the use of flavins and flavoproteins for novel photocatalytic, sensing, and biotechnological applications, 31,109−115 we expect that computational studies on these systems will continue to elucidate and guide experimental studies.In this work, we present a comparison of how different electronic structure methods describe the energetics of the low-lying excited states of flavin.
In the first part of this study, we compared S 1 excitation energies calculated with MR-PT2 methods and those computed with TD-B3LYP for flavin in different redox and protonation states.We focus on the S 1 state because, according to Kasha's rule, this is the state from which a system's photophysics and photochemistry typically occurs.An initial active space for CASSCF and MR-PT2 calculations is selected using autoCAS, but then the number of orbitals in the active space are gradually reduced to check the effect on the S 0 −S 1 excitation energy.Here, TD-B3LYP serves as a reasonable reference since it has been well tested against experimental spectra and shown to agree well when properly accounting for Franck−Condon factors. 21,38,39,61,67For most redox states, MS-CASPT2 (with the IPEA 0.25 shift) agrees well with TD-B3LYP excitation energies if a sufficiently large basis set is used.For the oxidized and semiquinone forms of flavin, which are planar, the convergence of the excitation energy is reached quickly with a small or moderate active space size.In the case of the hydroquinone forms of flavin, which are bent, a larger active space is needed for MR-PT2 methods to converge.For SS-CASPT2 and MC-PDFT, convergence with increasing active space size is not achieved for FlH 2 even when using a 16 electron and 14 orbital active space.
In the second part of this work, we compare TD-DFT, EOM-EE-CCSD, SOS-CIS(D), MC-PDFT, and MR-PT2 energy profiles along a PES connecting the S 0 and S 1 equilibrium geometries (i.e., tracing the Franck−Condon active mode) and along a PES connecting the S 1 (π,π*) and S 2 (n,π*) minima.Along the S 0 −S 1 path, we find that (TD)-B3LYP is in reasonably good agreement with MS-CASPT2 (IPEA = 0.25) energies, which explains its success in simulating UV/visible spectra.However, the same is not true for the S 1 −S 2 path, where TD-B3LYP has a significantly lower energy for the S 2 and T 2 states compared to MS-CASPT2 (IPEA = 0.25).TD-camB3LYP and TD-ωB97X-D, which are hybrid functionals that include a long-range correction that is missing in TD-B3LYP, appear to also have a reduced agreement with multireference methods.
While MS-CASPT2 with a minimal [2,2] active space and/ or 2-root state averaging can give VEEs that are in good agreement with larger active spaces, it is not recommended for studying regions far from the Franck−Condon point.
For all methods tested, the S 2 minimum lies above the S 1 minimum, suggesting that it is unlikely to be involved in flavin's photophysics.The T 2 (n,π*) state, on the other hand, may be populated by ISC from the S 1 (π,π*) state.The relative energies of these states are sensitive to the choice of computational method, and we posit that additional benchmarking is needed for a more quantitative description of ISC in flavins.Such benchmarks should ideally also test the effect of the method on the optimized geometry of the singlet and triplet (π,π*) and (n,π*) states.
−123 Some of this work is reviewed in ref 124.There, energies and wave functions of a wide range of methods were compared along paths relevant to PSB3's isomerization coordinate and S 1 excitedstate dynamics.Here, we again emphasize the usefulness of moving beyond vertical energies or energy differences between two points in benchmark studies.Instead, by comparing methods along a reaction or photophysical pathway, important differences in how these methods treat different regions of a potential energy surface may emerge.For PSB3, methods that could accurately reproduce the VEE could not necessarily describe the twisted PSB3 intermediates, which either displayed strong charge transfer character or diradical character.On the other hand, a few methods could treat both the planar and fully twisted PSB3 intermediates well but deteriorated when describing intermediate structures connecting the planar and 90°twisted geometries.Ultimately, through extensive benchmarks, methods were found that have a suitable error cancellation to treat different regions of PSB3's potential The Journal of Physical Chemistry B energy surface on a similar footing.This has enabled, for instance, the construction of model potentials suitable for running quantum dynamics for PSB3. 125Similarly, additional benchmark studies on flavin can lead to the adoption of methods that may be used routinely for studying the photophysics and photochemistry of this system in different protein environments.

Data Availability Statement
The main data from this work are presented in the manuscript and in the Supporting Information.Requests for additional data, including computational models and raw data files, should be directed to S.G. (sgozem@gsu.edu).
Cartesian (xyz) coordinate files of LF optimized in different redox and protonation states at the B3LYP/cc-pVTZ and MP2/cc-pVTZ levels of theory without symmetry constraints.Cartesian coordinates for structures along the S 0 −S 1 path and the S 1 −S

Figure 1 .
Figure 1.Structure of flavins in different redox and protonation states.The atom numbering is shown for the quinone (Fl) structure.

Figure 2 .
Figure 2. (TD)-B3LYP S 0 (brown), S 1 (blue), and S 2 (red) energies along two paths used for benchmarking electronic structure methods.The S 1 and S 2 labels are conserved throughout the PESs, despite state crossings, to reflect their energy ordering at the ground state equilibrium geometry.The larger circles indicate the minima optimized at the (TD)-B3LYP level of theory.Smaller circles are obtained by linear interpolation and extrapolation.At the bottom of each plot, a figure is shown overlaying the structures of all interpolated and extrapolated geometries.Bond thickness in this figure therefore correlates to the extent of change in coordinates.The labels above and below the structure indicate the bond lengths at the S 0 and S 1 (left) and S 1 and S 2 (right) structures, respectively.

Figure 3 .
Figure 3.Comparison of geometries and TD-B3LYP/cc-pVTZ excitation energies of flavin in five different redox states optimized with B3LYP and MP2.The schematic at the bottom left panel indicates the definition of the parameter d, which is difference in the maximum and minimum positions of heavy atoms along the z axis when the molecules are aligned as closely as possible along the x−y plane.The other five panels show the overlay of the structure for each redox state optimized at the B3LYP (blue) and MP2 (red) levels of theory.

Figure 4 .
Figure 4. CASSCF-DMRG[38,36] orbital entanglement diagram for oxidized LF.Gray circle size represents the degree of orbital entropy while line thickness indicates mutual information between orbital pairs.Red boxes highlight orbitals with high entropy.

Figure 5 .
Figure 5. 3-root MR-PT2 and 3-root state averaged CASSCF S 1 excitation energies as a function of active space [number of electrons, number of orbitals] for LF in the oxidized form.A 0.25 IPEA shift is used for MR-PT2 methods.The TD-B3LYP excitation energy is shown as a gray line for reference.XDW-CASPT2 energies are identical to MS-CASPT2 energies in this figure.

Figure 6 .
Figure 6.Orbital entanglement diagram for the semiquinone and hydroquinone reduced forms of LF.CASSCF-DMRG[39,36] actives spaces were used for the semiquinones and[40,36] for the hydroquinones.Gray circle size represents the degree of orbital entropy while line thickness indicates mutual information between orbital pairs.Red boxes highlight orbitals with high entropy.

Figure 7 .
Figure 7. 3-Root MR-PT2 and 3-root state averaged CASSCF S 1 excitation energies as a function of active space [number of electrons, number of orbitals] for LF in different redox states.The IPEA = 0.25 shift is used for all MR-PT2 methods.The redox state is indicated above each panel.The TD-B3LYP excitation energy is shown as a gray line for reference.XDW-CASPT2 energies are identical to MS-CASPT2 energies in these figures.

Figure 9 .
Figure 9. S 0 (bottom) and S 1 (top) energies along the S 0 −S 1 path reported relative to the corresponding S 0 minimum energy point along the path.Energies are shown for 3-root state averaged CASSCF, MR-PT2, MC-PDFT, TD-DFT, EOM-EE-CCSD, and SOS-CIS(D) methods, using colors indicated in the legend on the right.CASSCF and MR-PT2 energies employ a [14,12] active space.In the bottom panel, the circles indicate the interpolated geometry that has the lowest energy for each method.In the top panel, circles indicate the corresponding VEE (left side) and the AEE (right side) along this path.The latter is found as the point with the lowest S 1 energy along the path.XDW-CASPT2 energies are identical to MS-CASPT2 energies in this figure.
2 and T 2 (n,π*) states, give crossing geometries that are closer to the S 1 minimum, and have considerably smaller ΔE Sd 1 /Sd 2 and ΔE Sd 1 /Td 2 energies.EOM-EE-CCSD gives crossing geometries and ΔE Sd 1 /Sd 2 and ΔE Sd 1 /Td 2 energies that are intermediate between those of TD-DFT and MR-PT2.SOS-CIS(D), on the other hand, appears to significantly overestimate ΔE Sd 1 /Sd 2 and ΔE Sd 1 /Td 2 compared to the other methods.We note that EOM-EE-CCSD and SOS-CIS(D) calculations were carried out with a smaller basis set (double-ζ instead of triple-ζ), which may contribute to the difference observed.

Figure 10 .
Figure 10.PESs along the S 1 −S 2 path.In the left panel, we show the S 0 (bottom left panel) as well as S 1 and S 2 (top left panel) potential energies reported relative to the corresponding S 0 minimum energy point along this path.In the panel on the right, we show the S 1 and T 2 energies along the same path.Energies are shown for 3-root state averaged CASSCF, MR-PT2, MC-PDFT, TD-DFT, EOM-EE-CCSD, and SOS-CIS(D) methods, using colors indicated in the legend on the bottom right.MR-PT2 and CASSCF energies employ a [14,12] active space.In the bottom left panel, the circles indicate the interpolated geometry that has the lowest energy.In the tops panels, circles indicate crossing points between the S 1 state and the S 2 (left) or T 2 (right) states.XDW-CASPT2 energies are identical to MS-CASPT2 energies in these figures.

a
Those barriers are estimated for different electronic structure methods along the S 1 −S 2 path.The Journal of Physical Chemistry Bthe zeroth-order CASSCF wave function gives comparable S 0 and S 1 energy profiles along this path.However, using 2-root state averaging gives a qualitatively different S 1 energy profile.Similarly, using a minimal[2,2] active space, while successful at reproducing the VEE of flavin, has distorted S 0 and S 1 PESs when moving away from the S 0 equilibrium structure.This may result in artifacts far from the Franck−Condon region.■CONCLUSIONSMultiple computational studies have investigated the excitedstate electronic structure and photophysics of flavin.A nonexhaustive list includesrefs 21, 36−39, 50−81.

Figure 11 .
Figure 11.S 0 (bottom panel) and S 1 (top panel) energies along the S 0 −S 1 path reported relative to the corresponding S 0 minimum energy point along the path.Energies are shown for MS-CASPT2 with IPEA = 0.25 but with different number of roots used in the state averaging of the underlying CASSCF wave function.We also include one calculation where the minimal [2,2] active space is used instead of [10,10].

Table 1 .
NPEs for Different Methods Calculated Relative to the Reference XDW-CASPT2 PES

Table 2 .
Energy Barrier to the S 1 −S 2 Crossing Point (ΔE Sd 1 /Sd 2 ) and the S 1 /T 2 Crossing Point (ΔE Sd 1 /Td 2 ) Relative to the S 1 Minimum Energy a 2 path (21 structures each) (ZIP) Figures of the molecular orbitals included in the active space of CASSCF and MR-PT2 calculations; plots of CASSCF S 1 excitation energies as a function of active space and number of roots used in the state averaging for LF in the oxidized form (PDF)