Molecular dynamics (MD) simulations of System I. After proper setup, the system was then equilibrated in the following manner: (a) energy minimized using a combination of the steepest descent and conjugate-gradient methods; (b) heated up from 0 K to 300 K in NVT ensemble for 300 ps with a weak constraint of 25 kcal•mol -1 •Å -2 on the protein backbone atoms; (c) density equilibrated for 1 ns to achieve a uniform density at a target temperature of 300 K and a target pressure of 1.0 atm using an isothermal-isobaric Langevin thermostat and Berendsen barostat with a collision frequency of 2 ps -1 and a pressure relaxation time of 1 ps; (d) further equilibrated for 3 ns under the NPT ensemble with restraints removed to obtain a well-settled pressure and temperature. Finally, a 100 ns production MD simulation was performed. Nonbonded interactions were treated with Particle Mesh Ewald 73 method with a cutoff of 12 Å. Covalent bond containing hydrogen atoms were constrained using the SHAKE algorithm to enable an integration step of 2 fs. All MD simulations were performed with the GPU version of the Amber 22 package. Quantum mechanics and molecular mechanics (QM/MM) calculations. The final equilibrated geometry from the production MD trajectory was used as the initial input structure for the QM/MM calculations. QM/MM calculations were performed using ChemShell software, combining ORCA for the QM region and DL_POLY for the MM region with the AMBER force field. The QM/MM system contains the whole protein, counterions, and solvation water molecules within 8 Å of the protein. To include the polarizing effect of the MM region on the QM region, the electronic embedding scheme was employed in the QM/MM calculations. The time-dependent DFT (TD-DFT) method was employed for the QM region in excited state calculations. According to previous studies , the frontier molecular orbitals of HFOH mainly localize in the isoalloxazine ring and are not affected by the ribose phosphate side-chain. Therefore, only the lumiflavin group was included in the QM region to reduce computational cost without loss of accuracy, while other atoms were placed in the MM region. A hydrogen link atom with the charge-shift model was applied to treat the QM/MM boundary. Considering the charge transfer character of the first excited S 1 state and the well reproduction of the experimental fluorescence emission maximum (refer to Table ), the Coulomb-attenuated hybrid exchange-correlation functional CAM-B3LYP , which combines the features of the B3LYP functional with long-rang corrections using Hartree-Fock exchange, was used for the excited state calculations of the QM region. Geometry optimization of HFOH in its first excited S 1 state were calculated using a small basis set of 6-31+G** (labeled B1), with the residues within 6 Å of the chromophore in the MM region allowed to relax during the optimization process. The fluorescence spectroscopic properties were calculated with a larger basis set of def2-TZVP (labeled B2). Besides the method of CAM-B3LYP/def2-TZVP, four other functionals (B3LYP, M06-2X, 88 ωB97X-D3, and PBE0 ) and three other basis sets (6-31+G**, 6-311+G**, and 6-311++G**) were tested, which were summarized in Table . The fluorescence wavelength greatly depends on the choice of the functional, yet it is not sensitive to the choice of basis set. QM calculations. To ensure computational accuracy and efficiency, a simplified model of HFOH (4a-hydroxy-5-hydro-lumiflavin, denoted as LFOH) was used in QM calculations, where the ribose phosphate group was replaced with a methyl group (Scheme 2). The geometry of bare LFOH in its first excited state, used in System II, was optimized under vacuum using the TD-DFT method with the range separated CAM-B3LYP functional and the 6-31+G** basis set, at the same level for the QM region used in QM/MM calculations. The geometry used in System III was the QM/MM-optimized bare LFOH structure that accounted for the steric effects of the protein surroundings. The spectroscopic properties of LFOH, including emission energy E F , oscillator strength f, permanent dipole moment (μ) of the S 0 and S 1 states, and transition dipole moment μ 01 , were determined at the TD-CAM-B3LYP/def2-TZVP level, which is same with that used for QM region in QM/MM calculations. Vertical excitation energies were considered, ignoring vibronic effects. The effect of solvent polarity was investigated with System III using the polarizable continuum model (PCM) with five solvents (toluene, chlorobenzene, ethanol, DMSO, and water) covering a range of dielectric constant (ε). The effects of oriented external electric fields (OEEFs) were studied with System III using the planar QM/MM-optimized bare LFOH structure. The OEEFs were applied using the "Field = M ± N" keyword in Gaussian 16, which defines the axis of the electric fields, its direction along that axis (M) and its magnitude (N a.u., 1 a.u. = 5140 MV/cm). The short axis of isoalloxazine ring pointing from N 10 to N 5 was defined as the Z-axis and the long axis pointing from N 10 to N 1 was defined as the X-axis. (refer to Scheme 2). A range of moderate electric field strengths (F) between -0.0075 a.u. and 0.0075 a.u. were explored. A strong OEEF may change the energy order of the frontier orbitals of LFOH, which complicates the elucidation of the effect of OEEFs. It is noteworthy that the positive direction of the electric field vector is defined from the negative to the positive charge in Gaussian 16, which is opposite to the conventional definition in physics. To be consistent with the common sense of physics, the direction of the external electric field has been inverted in related discussions. The discussion of the effect of OEEFs are based on the results obtained using the TD-DFT/def2-TZVP method with the CAM-B3LYP functional.