Spectroscopic Properties of Lumiflavin: A Quantum Chemical Study

In this work, the electronic structure and spectroscopic properties of lumiflavin are calculated using various quantum chemical methods. The excitation energies for ten singlet and triplet states as well as the analysis of the electron density difference are assessed using various wave function‐based methods and density functionals. The relative order of singlet and triplet excited states is established on the basis of the coupled cluster method CC2. We find that at least seven singlet excited states are required to assign all peaks in the UV/Vis spectrum. In addition, we have studied the solvatochromic effect on the excitation energies and found differential effects except for the first bright excited state. Vibrational frequencies as well as IR, Raman and resonance Raman intensities are simulated and compared to their experimental counterparts. We have assigned peaks, assessed the effect of anharmonicity, and confirmed the previous assignments in case of the most intense transitions. Finally, we have studied the NMR shieldings and established the effect of the solvent polarity. The present study provides data for lumiflavin in the gas phase and in implicit solvent model that can be used as a reference for the protein‐embedded flavin simulations and assignment of experimental spectra.


INTRODUCTION
Flavins belong to a group of natural pigments responsible for the yellow color of various plants (flavus means yellow in latin). In 1937, Paul Karrer received the Nobel Prize in Chemistry for his characterization and synthesis of vitamin B 2 (riboflavin) (1). The yellow color of flavin in its oxidized form, as discussed in the present work, relies on the absorption with a maximum at 447 nm. A second electronic transition is in UV region and centered at k max of 360 nm (2).
The simplest flavin compound is the lumiflavin as shown in the oxidized and reduced forms in Scheme 1A,B. Substitution of a methyl group with an aliphatic side chain of 2,3,4,5tetrahydroxypentyl leads to riboflavin (Scheme 1C). In nature, enzymatic phosphorylation of riboflavin produces flavin mononucleotide (FMN), and further adenylation transforms FMN to flavin adenine dinucleotide (FAD) (3). These substituents do not have a significant effect on the absorption wavelength because they are not part of the conjugated system (2,(4)(5)(6), and thus, the present study is focused on the nonsubstituted lumiflavin.
The oxidized form of lumiflavin can be reduced twice. In the first reduction step, combined with a proton uptake, the neutral radical, also called semiquinone, is produced. Spectroscopically, the flavin radical can be distinguished from the oxidized form since it has absorption bands at 620, 580, and 500 nm, depending on the chromophore environment (5,7,8). Further reduction and protonation of the second nitrogen of the isoalloxazine ring leads to the form which has a major absorption at 250 nm and shoulders at 400 and 280 nm (9). Hence, flavins have the potential for transfer of electrons and protons, making them versatile for substrate modification and photochemically triggered redox reactions.
This abundance in reactivity makes flavins widespread functional cofactors in proteins. One class of flavoproteins is cryptochrome (10) closely related to DNA photolyases which conduct repair of the ultraviolet-damaged DNA (11,12). Both cryptochromes and DNA photolyases contain FAD as a cofactor and a conserved "tryptophan-triad" next to it. The latter forms an electron transfer chain when FAD is electronically excited and abstracts electrons to form an anionic radical (FAD •À ) (13). Also, BLUF domains, found in bacteria, contain FAD (BLUF = "bluelight using FAD") (14)(15)(16). They are involved in blue-light sensing. Phototropin is an FMN-containing protein controlling phototropism, that is, the orientation of plants toward light sources (17). The photosensory domain of phototropin includes two similar light-oxygen-voltage (LOV) domains, called LOV1 and LOV2. Each LOV domain contains an FMN chromophore within the protein binding pocket, bound through a network of hydrogen bonds and hydrophobic interactions (18). Obviously, the rich photo-and redox chemistry of flavins opens plenty of possibilities for biotechnological applications.
In order to interpret the experimental measurements and gain insight into the photochemistry of flavins, the use of the modern computational methods, in particular quantum chemical calculations, has proved to be vital. For example, previous theoretical studies by Zheng et al. (19) and Li et al. (20) have shown that the oxidized flavin adopts a planar structure whereas the reduced-form is bent. The calculation of vertical excitation energies has been helpful to assign the peaks in the absorption spectrum (21)(22)(23)(24)(25)(26). Also, the calculation of the vibrational frequencies has been instrumental in the assignment of the infrared spectrum (27,28).
In the present study, we re-examine the quantum chemical calculations of various spectroscopic properties. In the case of the vertical excitation energies of lumiflavin, a few lowerlying electronic states are reported often with contradicting energetic order or focus only on the bright excited states. Note that during the photochemical reaction, states which are dark initially might be populated as demonstrated recently by Falahati et al. (29). Further, a nonradiative decay from a singlet to a triplet state can occur via an intersystem crossing (30). The energetic order of the latter might be affected in TD-DFT calculations due to the triplet instability (31,32). Therefore, we would like to evaluate the energetic order of singlet and triplet excited states. Here, we report the results of calculations of 10 low-lying singlet and triplet excited states of lumiflavin using RI-ADC(2) and RI-CC2 as well as range-separated hybrid CAM-B3LYP methods with cc-pVTZ basis set. This work is focused on the analysis of the orbital character, charge-transfer, and excitation energies in both gas phase and solution.
Another focus of this work is the assignment of vibrational spectra. Previous studies have put an emphasis on the assignment of the carbonyl stretching modes. A large deviation between the experiment and quantum chemical calculations was described by Wolf et al. for these particular type of vibrations (27,28). The discrepancy was rationalized on the basis of environmental effects. However, in the present study we will evaluate the influence of anharmonicity. In addition, Raman and NMR spectra of the isolated and implicitly solvated lumiflavin are computed. The overall objective here is to review previous computational studies, revise some of the shortcomings and therefore create a reference for the gas phase that can be used in future studies where flavin is embedded within an environment.

MATERIALS AND METHODS
Geometry optimization. The geometry optimization for lumiflavin was carried out using the B3LYP functional (33) in conjunction with the 6-31G* basis set. The chemical structure of the above-mentioned flavin chromophore is shown in Figure 1A. All calculations were performed without a symmetry constraint. For DFT calculations, the D3 Grimme's dispersion correction was used (34). The B3LYP equilibrium geometry was used as our starting point for all calculations. The geometry optimization was carried out using in Q-Chem 4.4 program package (35).
Vertical excitation energies. The vertical excitation energies (VETs) of the ten low-lying singlet and triplet states were assessed using various wave function based methods as well as time-dependent density functional theory (TD-DFT). For all methods, the cc-pVTZ basis set and wherever required cc-pVTZ-RI auxiliary basis set were used.
Among the wave function methods, we have employed ADC(2) and CC2. The second-order algebraic-diagrammatic construction (ADC(2)) scheme (36) has favorable computational cost and good quality description of excited states. It has also the advantage of being size-consistent and having the eigenvalue equation formulated as a Hermitian matrix that is suitable for derivation of properties (37). The second-order approximate coupled cluster singles and doubles (CC2) method (38) is also proved to be efficient for the accurate description of excited states. The accuracy of CC2 was assessed previously in the benchmark studies by Schreiber et al. (39) using 28 organic molecules. Recently, also an assessment of the ADC(2) method was made by Harbach et al. (40) using the same set of molecules. The singlet excitation energies of CC2 method were found to have an accuracy of 0.29 AE 0.28 eV (mean value AE SD) for 103 vertical excitation energies compared to the theoretical best estimates, while ADC(2) produced 0.22 AE 0.38 eV for 104 vertical excitation energies. Similarly, for triplets a comparison of 63 vertical excitation energies showed a deviation of 0.17 AE 0.13 eV and 0.12 AE 0.16 eV for CC2 and ADC(2) methods, respectively. The two wave function-based methods were used in the "resolution of the identity" (RI) formulation (41), which is an order of magnitude more efficient in comparison to the conventional algorithms. The calculations for RI-ADC(2) and RI-CC2 were performed in Turbomole 7.0 (42).
Time-dependent density functional theory TD-DFT (43) method is employed for calculation of electronically singlet and triplet excitation energies. Two functionals were involved: B3LYP and the long-range corrected version: coulomb-attenuating method (CAM)-B3LYP (33,44). The hybrid functional B3LYP is known to give accurate nÀp* and pÀp* excitation energies, with an expected error range of ca. 0.20-0.25 eV from the experimental values (43). A benchmark study of TD-DFT by Silva-Junior et al. (45) revealed a deviation of À0.07 AE 0.33 eV for singlet and À0.45 AE 0.49 eV for triplet vertical excitation energies, with reference to the theoretical best estimates reported in (39). Notably, these statistics are based on consideration of 104 singlet and 63 triplet excitation energies, which further infers that the hybrid TD-B3LYP describes the excitation energies more accurately compared to other functionals like BP86 and BHLYP (45). Further, TD-CAM-B3LYP method can correctly describe excited states with charge transfer (CT) character (46). Another important issue is the "triplet instability problem" in the ground state wave function, as described by Peach et al. (31). In their work, they have shown that the use of hybrid and range-separated functionals within the "full" TD-DFT scheme does not produce accurate singlet-triplet energy gaps. Instead, the use of Tamm-Dancoff approximation (TDA) was reported to show a dramatic improvement of the singlet-triplet gaps (31,32). Hence, in order to have more accurate results, we have used TDA for the calculation of IR and Resonance Raman spectra. One of the most valuable tools in the characterization of the molecular structure is the vibrational spectroscopy, which has two main directions: Infrared and Raman. The main difference between the two types lies in the selection rules of the transitions. The vibrational modes which change the dipole moment of the molecule are active in the IR, whereas the modes which change the polarizabilities are Raman active. In contrast to the conventional Raman spectroscopy, which is based on the direct absorption of the electromagnetic radiation in the IR region, the resonance Raman (RR) method is based on the vibronic transitions (which occur between different electronic and vibrational states) that can lead to a significant enhancement of the Raman bands (50). From the computational perspective, the vibrational modes correspond to energy levels that are obtained from the solution of the nuclear wave function equations. Hence, the fundamentals for the IR and Raman signals are computed in the same way but the calculation of the intensities is different. The IR intensities require mixed second derivatives of the energy with respect to geometric displacement and an external electric field, while the Raman intensities require mixed third derivatives to account for a change in the molecular polarizability along the geometric motion.
The gas-phase IR intensities were computed using the B3LYP functional and the 6-311G* basis set in Gaussian 09 program package. The obtained harmonic frequencies were scaled by a factor of 0.966 as suggested in the Computational Chemistry Comparison and Benchmark Database (51). In addition, we have tested the effect of the anharmonicity by using the Generalized second-order Vibrational Perturbation Theory (GVPT2) (52).
The Raman spectrum was calculated using the ORCA program at the B3LYP/6-311G* level of theory. For the RR spectrum, the gradients were obtained by the TD-B3LYP method. From previous benchmarking studies (53), B3LYP functional is found to be the most efficient for determination of the ground state properties (vibrational frequencies and normal modes) as well as for the resonance Raman (RR) intensities. The method employed in this work is based on Independent mode, displaced harmonic oscillator (IMDHO) as described by Petrenko and Neese (54). The solvation effect was included in the calculation of the RR spectrum, using CPCM. The ORCA modules "orca_mapspc" and "orca_asa" were used for extracting the polarizabilities and simulating the RR spectrum, respectively (48). For the visualization, Lorentzian function with Full width at half maximum (FWHM) of 10 cm À1 is used for the peak broadening.
NMR spectrum. The chemical shielding values were calculated using GIAO approach in Gaussian 09. These calculations were performed at B3LYP level using the 6-311G* basis set. The polarity effects of water and cyclohexane were taken into account using CPCM.
Experimental section. For measuring the UV-vis absorption spectrum, lumiflavin was dissolved in chloroform and the spectrum was recorded from 250-600 nm by UV-2450 spectrometer (Shimadzu). The data acquisition is performed with a step interval of 0.5 nm.
The infrared absorption spectrum was measured in this work. For this purpose lumiflavin powder was first mixed with KBr powder (wL/ wK = 1:50) and grinded until they became homogeneous. Subsequently, the mixture was subjected to high pressure for two minutes to form a thin transparent pellet. The pellet was then measured by Fourier transform IR spectrometer (Thermo Nicolet) from 4000-400 cm À1 at a spectral resolution of 1 cm À1 . The data processing was done using Omnic32 software.

Optimized geometry and IR spectrum
The optimized geometry of the oxidized form of lumiflavin along with the bond lengths is given in Figure 1A. The numbering of atoms is used consistently throughout the manuscript. The vibrational frequencies for bond stretching modes are related to the mass of the atoms and to the strength of the corresponding chemical bond involved. The IR intensities were calculated in vacuum at the B3LYP/6-311G* level of theory. A scaling factor of 0.966 was applied to the harmonic frequencies according to the employed methodology (Table 1) (51). To the best of our knowledge, there is no experimental gas-phase IR spectrum of lumiflavin. Hence, we have measured a spectrum in the KBr disk and compared the observed peaks with the calculated frequencies. The solid environment present in the crystalline KBr pellet provides sharp peaks, but the tight packing in the solid could lead to interactions with lumiflavin and in addition the thickness of the pellet could affect the absorbance of the sample. Therefore, we cannot expect a direct correspondence with the calculated frequencies of the isolated lumiflavin even after the scaling of the frequencies. The assignments of the peaks in Figure 1B are made using a reference study based on isotope labeling (55). In the following discussion, we focus mainly on the most intense bands of the spectrum.
Two amide I vibrations (C 4 =O 4 and C 2 =O 2 stretching modes coupled with N 3 -H bending) in the experimental solid-state spectra are observed at 1708 and 1663 cm À1 . Similar values were obtained for the flavin in the protein environment: for FMN (in dark state) in LOV2 domain at 1714 and 1678 cm À1 , as determined experimentally from difference spectrum (56) and in BLUF domain at 1710 and 1641 cm À1 , based on difference spectrum (57) and isotope labeling (58). Further, for the oxidized lumiflavin in LOV1 domain (C57S mutant) these vibrations were assigned to 1693 and 1681 cm À1 , as determined from the difference spectrum from vibrations of dark and light state (59). The calculated vibrational modes of the C 4 =O 4 and C 2 =O 2 stretching appear at 1738 and 1727 cm À1 , respectively. Although a scaling factor of 0.966 has been applied the deviation, especially for the C 2 =O 2 mode, is fairly large. However, it is not surprising that one global scaling factor is not sufficient to have a close agreement with the experiment. Mroginski et al. (60) have reported the use of 11 scaling factors for the linear tetrapyrrole chromophore in phytochromes in order to achieve an agreement with the experiment of ca 10 cm À1 . In the specific case of the two C=O stretching modes in flavins, Wolf et al. (27) have first proposed hydrogen bonding to be responsible for the discrepancy between the experiment and the quantum chemical calculation. In a subsequent study, the same group has attributed this discrepancy to the bulk solvent effect (28). By adding an implicit solvent to the calculation the agreement with the experiment in DMSO has improved for the two C=O stretching modes. Nevertheless, the measured C=O stretching frequencies in DMSO are in good agreement with our experimental results obtained in the solid state: 1711 and 1676 cm À1 . Hence, this agreement might question the interpretation of the solvent. The deviation could be also due to the lack of anharmonicity in the calculated frequencies. In order to test this hypothesis, we accounted for the missing anharmonic correction by the GVPT2 method (see Supporting Information). The obtained anharmonic frequencies are found to be in a better agreement with the experimental counterparts for the C 4 =O 4 and C 2 =O 2 vibrations (1709 and 1695 cm À1 , respectively). However, the GVPT2 method has downshifted all frequencies which were previously in good agreement with experiment. Thus, our findings point to a possible environmental effect which is the major reason for the difference observed in the C=O stretching mode frequency as reported previously by Wolf et al. (28). However, the environmental effect is not specific to the polar solvent DMSO.
The contributions from the skeletal stretching of ring I coupled with ring II and III in the solid-state spectrum were found at 1621 and 1461 cm À1 . These vibrations appear at 1615 and 1421 cm À1 in the computed spectrum. Further intense bands in the measured spectrum are found at 1583, 1552 and 1169 cm À1 , which also correspond to ring II stretching vibrations (55). The gas-phase calculation yields bands at 1570 and 1535 cm À1 that correspond to the ring II stretching with main contribution from N 10 -C 9a and C 4a -N 5 stretching. These ring II stretching modes were also described in the experimental IR difference spectrum of flavin from plant cryptochrome (CRY1) at 1578 and 1543 cm À1 (61). Also, in earlier reports (59,62), similar vibrations, originating from ring II of the isoalloxazine were assigned to 1550 and 1584 cm À1 based on analysis of difference spectrum. The calculated vibrational band at 1151 cm À1 corresponds to the C 4a -C 4 and N 5 -C 5a stretching coupled with N 3 -H bending. A low-intensity N 3 -H bending mode measured at 1413 cm À1 was computed at 1375 cm À1 . The large deviation can be explained by the involvement of the N 3 -H group in hydrogen bonding within the protein binding pocket (63,64). For example, in FMN bound to LOV domain, the N 3 -H bending mode appears at 1443 cm À1 as obtained from quantum chemical calculations (65). We also noted a delocalized skeletal vibration accompanied with N 10 -C 9a and C 4a -N 5 stretching in our calculated spectrum at 1284 cm À1 . These vibrations were also referred to ring II stretching modes, which correspond to the peak at 1301 cm À1 in the experimental solid-state IR spectrum of lumiflavin. However, according to our calculations this also involve ring III. The amide III vibration (C 2 -N 3 antisymmetric stretching mode) coupled with N 3 -H bending mode was found experimentally at 1238 cm À1 , which in the gas phase occurs at 1169 cm À1 . A similar antisymmetric stretching is assigned experimentally through difference spectrum in FMN bound to LOV2 domain at 1395 cm À1 , which also involves C 4 -N 3 vibration along with stretching vibrations from C 2 -N 3 and bending  (56)  1169 1151 C 4a -C 4 , N 5 -C 5a , N 3 -H bending vibrations from N 3 -H (56). Frequencies lower than 1000 cm À1 are mostly referred to the skeletal deformation vibrations (55,63). For example, a calculated peak at 800 cm À1 , corresponds to a stretching vibration of rings I and II. In this frequency region the vibrations are mostly delocalized over the fused multiple rings (66).

Raman and resonance Raman (RR) spectra
The computed and measured Raman as well as RR spectra of the oxidized lumiflavin is shown in Figure 2. The experimental Raman spectrum of lumiflavin was measured in powder (67), since we cannot account for this environment we have chosen to calculate the Raman spectrum in vacuo at the B3LYP/6-311G* level of theory (Figure 2A,B). However, we have selected the experimental RR spectrum obtained for lumiflavin in the riboflavin-binding protein (RBP) in water solution that was excited at 488 nm (68). Therefore, we decided to compute the spectrum in the implicit water model (CPCM) at the B3LYP/6-311G* level of theory ( Figure 2C,D). Hence, not only the intensities but also the frequencies will be different between the computed Raman and RR spectra. The frequencies of the calculated Raman spectrum of lumiflavin are in qualitative agreement with the experimental results in the region 1700-1100 cm À1 , given that the latter is measured in powder (67). The band at 1080 cm À1 in the experimental spectrum corresponds to stretching vibration of C 7 -C 7a and C 8 -C 8a along with skeletal deformation of rings I, II and III, which was found at 1094 cm À1 in the calculated counterpart. The experimental peaks within the fingerprint region of 1300-1100 cm À1 correspond to a combination of the bending of C 6 -H,  1368 cm À1 , respectively. Peaks in the experimental spectrum at 1415 and 1465 cm À1 are due to the ring II and III vibrations, C 7a and C 8a methyl deformation, as well as C 9 -H bending, which are shown at 1402 and 1469 cm À1 in Figure 2A. The largest discrepancy is found for the three intense peaks at 1574, 1588 and 1624 cm À1 in the calculated spectrum. Their relative intensity is much higher than in the experimental spectrum. They are best assigned to the peaks measured at 1552, 1583 and 1623 cm À1 . The first mode is a skeletal stretching of ring II and III with contributions from N 1 -C 10a and the C 11 -methyl group, the second is associated with the ring III stretching, C 7a and C 8a methyl groups and the third one is stretching between N 1 -C 10a -C 4a -N 5 , coupled with C 8 -C 9 . The stretching mode of C 4a -N 5 is also found at 1580 and 1585 cm À1 in case of FMN in LOV1 domain and in water (62).
As mentioned in the Methodology section, the intensity in the RR spectrum is greatly enhanced for those modes that are active in the excited state. We will focus on the most intense peaks in the following discussion because we have described the assignment of the Raman peaks already above. The calculated frequencies in the Raman and RR spectra should be identical and only the relative intensities are altered. However, we have used an implicit water solvent in the RR calculation because the experimental spectrum was also recorded in water solution of RBP. One cannot expect a quantitative agreement with the experiment since we do not have the specific protein interaction, nevertheless, the qualitative trend can be observed. In the computed spectra, we have noted an enhanced intensity when going from Raman to the RR spectrum at 1345 and 1368 cm À1 (Figure 2A), while the two peaks at 1217 and 1240 cm À1 have decreased intensity ( Figure 2C). The latter two peaks originate from ring I vibrations: the one at 1217 cm À1 due to C 2 -N 3 and C 4 -C 4a stretching and the other one at 1240 cm À1 due to the bending of C 2 -N 3 -C 4 . They agree well with the experimental counterpart, where the most intense peak is found at 1243 cm À1 ( Figure 2D). Further, intense RR peaks in calculated spectrum are found at 1541, 1604 and 1663 cm À1 with two of them corresponding to the intense peaks at 1574 and 1588 cm À1 in the computed Raman spectrum and are described above. Due to the lack of the baseline correction in the experimental counterpart, their assignment is difficult. In the measured RR spectrum the signal at 1553 cm À1 corresponds to the peak at 1541 cm À1 of calculated RR and is also found at 1553 cm À1 in a RR spectrum of FAD in water (69), while the peak at 1585 cm À1 corresponds to 1604 cm À1 in the computed RR spectrum. The intense peak at 1663 cm À1 of Figure 2C corresponds to a less prominent peak at 1670 cm À1 in the calculated Raman and to 1630 cm À1 in the experimental RR spectrum. It is associated with the stretching of C 10a -C 4a , C 9 -C 9a , C 5a -C 6 and C 6 -C 7 from ring II and III.
In the extended low-energy region of the RR spectrum (400-1100 cm À1 ), that is not covered in the Raman spectrum, the most intense peak is found at 523 cm À1 . It corresponds to 527 cm À1 in the experimental RR spectrum ( Figure 2D) and is attributed to the stretching of ring I and II as well as the C 7a and C 8a methyl groups (55,66). Other peaks at 604 and 640 cm À1 that agree well with the experimental RR spectra arise from the bending of the C 2 -N 3 -C 4 and N 5 -C 5a -C 6 angles, respectively. The computed RR peak at 748 cm À1 corresponds to the vibration of ring III, which is shifted in the experimental counterpart by 7 cm À1 . Further peaks at 793 and 839 cm À1 , correspond to the displacement of C 6 -C 5a along with stretching between C 2 -N 3 shifted by -3 cm À1 and C 4a -N 5 -C 5a shifted by 4 cm À1 , respectively. The vibrations at N 10 were assigned to the intense signal at 997 cm À1 in the experimental spectrum, which was found at 997 cm À1 in the calculated one.

Assessment of vertical excitation energies of lumiflavin
A comprehensive understanding of the electronic structure associated with the absorption peaks of the lumiflavin spectrum is the main target of numerous computational studies (21,22,(24)(25)(26)(70)(71)(72)(73). The nature of the excited states and their energetic order is also important for the study of the photochemical reaction initiated at the Frank-Condon point. Therefore, the aim of this work is to assess the excited states in lumiflavin using ab initio quantum chemical methods ADC(2) and CC2 and compare the results to previous studies.
We have calculated the 10 lowest singlet and triplet excited states which can become relevant in the investigation of the excited state reactivity including singlet-triplet transitions.
A summary of vertical excitation energies of lumiflavin is presented in Table 2. In the analysis of the singlet excited states, we consider a bright state to have an oscillator strength (f) <0.1. Molecular orbitals (MOs) that are involved in the low-lying electronic transitions are shown in Figure 3. The first bright state (S 1 ) obtained with CC2 and ADC(2) is in the visible spectral range (ca. 400 nm), whereas other bright states fall into the near-UV spectral region. In general, the vertical excitation energies computed with ADC(2) are typically lower by ca 0.10 eV for both singlet and triplet transitions, compared to that of CC2. In ADC(2), four bright states are found among the 10 lowest singlet states: 2.94 eV (S 1 ), 4.06 eV (S 4 ), 4.79 eV (S 7 ) and 4.91 eV (S 9 ). However, CC2 calculation reveals three bright states located at 3.06 eV (S 1 ), 4.13 eV (S 4 ) and 4.97 eV (S 9 ). This is because S 7 has an oscillator strength of 0.088 at the CC2, which is close to the arbitrary selected cutoff of 0.1 for the bright states. Therefore, the only change in the relative order of the first 10 singlet states is found to be S 7 at CC2 level of theory, which corresponds to S 9 at ADC(2) level of theory. We attribute this swap in the position to the close proximity of the states S 7 , S 8 and S 9 , which are found within 0.04 eV at CC2 level and 0.12 eV at the ADC(2) level, which are nearly degenerate. However, the nature of singlet excited states for CC2 and ADC (2)  . Another high-level method is XMCQDPT2; however, the published results include only two pÀp* states (25). While the DFT/MRCI energies are slightly lower than CC2, the XMCQDPT2 energies are slightly higher. However, the reported XMCQDPT2 values by Udvarhelyi et al. (25) were obtained for a geometry that was optimized at a different level of theory (Complete Active Space Self-Consistent Field). The authors of that study have reported a lowering of the excitation energies in case the geometry is optimized at the B3LYP level of theory resulting in close agreement with CC2 values in the present work ( Table 2).
The experimental results show that the first peak of lumiflavin in water appears at 2.79 eV and in the protein environment at 2.77 eV (71,74). A comparison of this peak with the first lowlying bright state of CC2 in gas phase yields a blueshift of 0.27 eV. The second experimental peak is found at 3.38 and 3.77 eV in water and protein, respectively, which in gas phase is blueshifted by 0.75 and 0.36 eV. The assignment of the third experimental peak with reference to the gas-phase results is complicated because several states with non-negligible oscillator strength are found in a close proximity. In case of S 6 and S 7 states, the oscillator strength is 0.016 and 0.088, respectively. However, the S 9 state is found at 4.97 eV that is within 0.04 eV  of S 7 and has a nearly 10-fold increased oscillator strength (f = 0.776). Hence, we expect the brightest peak to overlap with S 7 and impede the assignment of the experimental spectrum. Due to the high oscillator strength of S 9 and the high intensity of the third peak ( Figure 5), we have given the values in the same row of Table 2.
For a detailed characterization of the excited states, we have analyzed the electron density differences (EDD) between the S n and S 0 . The EDDs were calculated for the four excited states with the highest oscillator strength from CC2 calculations (Figure 4). The ground and excited state densities are delocalized over the entire molecule, except for the S 0 -S 4 transition, which possess a charge-transfer from ring III to rings I and II. The S 1 -S 0 EDD is more pronounced on the rings I and II, which shows that the electron density is transferred from N 1 to N 5 upon excitation. These results are in agreement with the previous reports (72).
In comparison to the results of the wave function methods, the previously reported TD-DFT results by Neiss et al.  Information Table S2), suggest not only the different ordering of the singlet but also triplet excited states. The TD-DFT calculations using the B3LYP functional result in a bright, first singlet state with pÀp* character. However, the following three states are dark nÀp* states in calculations with a double-zeta basis set. Only the fifth excited state corresponds to the bright pÀp* transition. If a triple-zeta basis set is employed, there are only two of these dark states between the two bright states in agreement with ADC(2), CC2, and DFT/MRCI. Also the long-range corrected functional CAM-B3LYP can give accurate description for state ordering among the low-lying excited states. This has been shown for double-zeta basis set by List et al. (73) as well as triple-zeta basis set in Table 2. The fact that CAM-B3LYP is adjusting the energetic order indicates a partial charge transfer character of the S 4 state which is not correctly described using B3LYP. The nature of the charge transfer is evident from difference density plots for the bright states in Figure 4, as discussed above. A detailed inspection of the excited states shows that despite of the same order of bright states as in CC2, the electronic character for two pairs of states is interchanged at the CAM-B3LYP level (S 2 and S 3 , S 7 and S 8 ). The states at 3.55 eV (S 1 ), 4.37 eV (S 4 ), and 5.42 (S 9 ) are found to be bright which is consistent with the CC2 results. The highest computed excited state S 10 involves different orbitals compared to CC2, and therefore, is omitted from Table 2. The range separation in CAM-B3LYP correctly describes the order of the states; however, the excitation energies are overestimated because the cparameter (controlling the short-and long-range separation) is not universal (76).
In the case of the triplet excited states, B3LYP yields four states (in work of Falkl€ of et al. (77) three states, but S 1 and T 4 are nearly degenerated) below the first singlet. For TD-CAM-B3LYP, there are three triplet states and for DFT/MRCI two triplet states below S 1 (Table 2). Both, ADC(2) and CC2, methods indicate that only one triplet state is below the lowest singlet excited state, which is in-line with CIS(D), a simple configuration interaction method (Table S2). Hence, the problem with triplet states could be attributed to the documented triplet instability in TD-DFT (31,32). The S 1 -T 1 energy gap is 0.54 and 0.50 eV for ADC(2) and CC2, respectively. In summary, we conclude that the use of TD-DFT should be avoided in particular for studies of the lumiflavin photochemistry involving singlet-triplet transitions.

Solvatochromic effects on vertical excitation energies
In the previous section, we have shown that the first two bright states from ADC(2), CC2, and CAM-B3LYP calculations involve p H ? p* L and p H-1 ? p* L transitions, and their positions are found in qualitative agreement with the experimental peaks (Table 2). However, the gas-phase results are blueshifted with respect to the absorption peaks of lumiflavin measured in water and within the protein embedded system (73). Our recorded spectrum in chloroform is shown in Fig (78,79). In this section, we want to address to the solvent effect on the few low-lying excited states that are needed to understand the experimental spectrum (80,81). In order to analyze the solvation effect in a systematic way, we calculated the VETs using the Conductorlike Polarizable Continuum Model (CPCM) (47), with water, chloroform, and cyclohexane solvents. Unfortunately, the CPCM implicit solvent model is not available for CC2 or ADC(2) methods in Turbomole. Thus, we have decided to use the CAM-B3LYP method to study the effect of solvation ( Table 3). Comparison of the low-lying bright states in the gas phase with those computed in water reveals a redshift from 0.14 to 0.28 eV. This solvatochromic shift improves the gasphase results in comparison to the experimental peak maxima (Tables 2 and 3). If we assume the shift is independent of the electronic structure method, then the solvent model will also improve the CC2 excitation energies. In case of cyclohexane, the shift of the bright state energies with respect to the gas phase is marginal (ca. 0.1 eV), due to the low dielectric constant of this nonpolar solvent.
Calculations in other solvents revealed different solvatochromic effect on the bright states. The series of the solvents includes cyclohexane, chloroform, and water. The polarity of the solvent is increasing from cyclohexane having the lowest value that makes it the closest to gas phase to water with highest polarity. The S 1 energy is nearly unaffected but the oscillator strength is decreasing when going from the nonpolar cyclohexane to the polar water. In this series of solvent, the S 4 shows a decrease in excitation energy and an increase in oscillator strength. The largest change in the excitation energy is found for the S 6 , where a blueshift of ca. 0.5 eV from cyclohexane to water is observed. In comparison to the experiment, the computed values (chloroform) are systematically blueshifted by ca. 0.65 eV. This might be due to simple nature of the implicit solvent model, which for example cannot account for hydrogen bonding. To illustrate, List et al. have performed a simulation of the absorption spectrum of lumiflavin with polarizable embedding technique to account for the solvent effects. The obtained absorption maxima were blueshifted by ca. 0.4 eV in water (73,74), and the order of states was slightly different: the dark S 2 state is destabilized and shifted above S 4 . The consideration of low energy, low-frequency modes is another important factor as demonstrated via sampling of the ground state PES by Marazzi et al. (82).
Vibrational broadening can be also included via the vertical and adiabatic Franck-Condon approaches (VFC and AFC) as demonstrated by Karasulu et al. for flavin derivatives (83). However, since the gas-phase experiment on lumiflavin has not been done yet, the reproduction of the bandshapes of UV spectra is out-of-scope of the current study. According to the reference above, inclusion of the vibrational effect will provide a redshift. We have observed that CAM-B3LYP excitation energies in PCM solvent overestimate the measured absorption maxima by 0.64 eV for S 1 . Hence, the vibrational effect will reduce the difference between the computed and measured spectra. However, in this study we have mainly focused on the relative position of the bands.
Note that the low-energy band consists of three peaks, which yet correspond to the transitions to the different vibrational levels of the same electronically excited state ( Figure 5). This was previously proven by Klaum€ unzer et al. (84), by simulating a vibronic spectrum of riboflavin. Also, it can be seen from our computed transitions: the energy gap between first two bright states is 0.74 eV (CAM-B3LYP, chloroform), which is significantly higher than the gap between different shoulders of the first absorption band (ca. 0.15 eV).

NMR chemical shielding analysis
The NMR chemical shielding values are given in Table 4. The experimental data for the carbon and nitrogen atoms were obtained from previous studies of the oxidized form of tetraacetylriboflavin (85,86). The conversion of the measured chemical shift values to chemical shieldings is described in ref (86). Since chloroform can only donate a hydrogen bond but not accept it, Walsh and Miller (86) have corrected the values for N 1 and N 5 by À12.5 ppm. Further, it is also important to note that the experimental shifts were measured in chloroform. Hence, in the present work we have investigated the solvent effect on the shielding values by choosing three difference solvents with increasing polarity. The three solvents are cyclohexane, chloroform, and water consistent with the excitation energy calculations.
The calculated NMR shielding values are sensitive to the choice of basis set. A large shift in these values is observed when going from double-zeta (6-31G*) to triple-zeta (6-311G*) basis set (Table 4 and Table S3). The aim of our study is not to reach a quantitative agreement but rather compare the relative  carbon and nitrogen shieldings which agree well with the experimental trends for both basis sets. Among all carbon atoms, the highest deshielded values belong to C 7a and C 8a . The downfield shift is attributed to the fact that the electron density of these methyl groups at ring III shifts toward rings I and II which contain electronegative atoms (87). The introduction of implicit solvent in the calculation does not change the values of these carbon atoms. A significant effect of the increased polarity when going from the gas phase to the water environment is found in the upfield shift of C 2 , C 4 , C 7 , and C 8 by 4.92, 4.7, 5.23 and 6.57 ppm, respectively. In the nonpolar solvent, the maximum upfield value is found for C 8 with 2.35 ppm. Among the nitrogen shielding values, only the one of N 5 is negative. A negative value indicates an electron-rich atom, which in case of N 5 is due to the lone pair and the strong shielding effect of the fused ring current. Also, N 5 is further away from the electron withdrawing carbonyl group compared to N 1 that also has a lone pair. The effect of polarity (water) is responsible for the highest deshielding noticed for N 1 with a downfield shift of 17.81 ppm, while the other three nitrogens are not affected. This change is less pronounced in the cyclohexane solution, namely by 6.79 ppm. This indicates a high sensitivity of N 1 to the environment. In contrast, the shielding at N 10 is upshifted. The change in water is by 13.44 ppm and in cyclohexane by 4.86 ppm. In summary, based on the assignment of the NMR shieldings we confirm the experimental findings that identified N 5 as the most favorable nucleophilic reaction center (87,88). However, we find the shielding of N 1 to be the most sensitive to the polarization of the solvent environment.

CONCLUSIONS
This work provides a detailed investigation of the various spectroscopic properties of lumiflavin in the gas phase and in solution to give a reference for further photochemical studies. Ten low-lying singlet and triplet states were calculated by various levels of theory. The nature of all excited states was analyzed in terms of the character of the involved orbitals and electron density difference for the bright states. Furthermore, the IR, Raman, as well as resonance Raman spectra were simulated in order to assign the experimentally observed vibrational spectra. In addition to the above, the GIAO NMR shielding of nitrogen and carbon atoms was computed and compared to the experiment.
The vertical excitation energies were computed using the accurate ab initio methods CC2 and ADC (2). Both of them yield the same order of the singlet excited states as previously reported DFT/MRCI, nevertheless, the ordering of the triplet states is qualitatively different. The results of CC2 and ADC (2) were also compared to the results of the earlier TD-DFT studies. It was found that the most popular B3LYP functional does not give an accurate description of the relative state ordering for singlet and triplet states. However, the inclusion of range separation (CAM-B3LYP) provides a more balanced description of the singlet excited states, yet, the order of the triplets is incorrect. Thus, on the basis of the current and previous studies, we make several conclusions: (1) Despite the identification of the two lowest bright states of lumiflavin, there are two dark nÀp* states inbetween which can cross in the course of the photoreaction and thereby affect the reaction; (2) TD-DFT methods provide reasonable excitation energies for the bright states; however, the relative ordering of the states and singlet to triplet transitions might not be correctly described. (3) As mentioned above, CAM-B3LYP is free of the limitations of B3LYP for the singlet states. Further, our study of the solvatochromic effect, revealed no dependence of the first bright excited state on the solvent polarity, whereas remaining bright states undergo a significant change in both excitation energies and oscillator strength. Despite the fact that deviation from the experimental values is significant (0.64 eV), it is systematic.
The comparison of the simulated and the measured IR spectrum showed good agreement for the majority of the frequencies. However, two modes, namely C 4 =O 4 and C 2 =O 2 show a significant deviation by ca. 30 and 64 cm À1 , respectively. The account for anharmonicity reduces one frequency more than the other, that is difference of 1 and 32 cm À1 . However, the agreement of the most of the other frequencies deteriorates. Hence, we conclude that the lack of anharmonicity is not the major source of the large deviation in the carbonyl stretchings because this correction shifts most of the vibrational frequencies. Further calculations included the conventional and resonance Raman spectra where we found good agreement to the corresponding experiments in solid state (Raman) and in protein (RR). In particular, the enhancement of the Raman intensity due to the resonance condition was qualitatively reproduced. However, the identification of the intensities in the region 1200-1700 cm À1 in the experimental RR spectrum was hindered by the lack of the baseline correction.
Finally, we studied the GIAO NMR shielding values of the carbon and nitrogen atoms in lumiflavin. We have identified which shielding values are sensitive to solvent polarity and managed to observe trends that help to assign experimental spectra. The values of carbon atoms C 2 , C 4 , C 7 , and C 8 changed the most while for the nitrogens it was N 1 and N 10 .
The importance of this work lies in the identification of the spectroscopic values which are difficult to observe (e.g. dark electronic states) or to assign unambiguously (e.g. solvatochromic effects). Overall, we believe that this computational work will serve as a reference for embedded calculations and will also support the assignment of experimental spectra of flavins.

SUPPORTING INFORMATION
Additional supporting information may be found online in the Supporting Information section at the end of the article: Table S1. Experimental vibrational frequencies compared with calculated anharmonic frequencies. Table S2. Vertical excitation energies computed with RI-CIS (D) and TD-DFT B3LYP. Table S3. Vertical excitation energies computed with TDA TD-DFT B3LYP with CPCM solvent. Table S4. GIAO NMR shielding values computed at B3LYP/ 6-31G* level.