Simulating accurate condensed phase spectral lineshapes requires the inclusion of both vibronic and solvent effects. Such effects contribute to spectral width, as well as shape, with vibronic transitions often causing a high energy tail in absorption spectra and solvent effects often leading to inhomogeneous broadening. Although this inhomogeneous broadening is generally accounted for with a phenomenological broadening parameter, it is specific interactions between a chromophore and the solvent, such as solvent hydrogen bonding and dipole-dipole interactions, that lead to the energy fluctuations responsible for inhomogeneous spectral broadening. Sampling from an ensemble of independent chromophore-solvent nuclear configurations can capture environmental effects within a nuclear ensemble approach to simulating the lineshape, but this nuclear ensemble approach lacks vibronic transitions. The vibronic effects on a spectrum can be accounted for with a Franck-Condon calculation, often performed within the harmonic approximation with an implicit polarizable continuum model for the environment. Alternative methods can simultaneously include vibronic and specific environmental effects. Recent work by some of the authors shows how to combine the Francka) These authors contributed equally to this work. b) Electronic mail: lshi4@ucmerced.edu c) Electronic mail: cisborn@ucmerced.edu Condon approach with an ensemble of explicit solvent configurations. An alternative approach is to follow the dynamics of the chromophore and solvent by saving a molecular dynamics (MD) trajectory of chromophoresolvent configurations for the computation of energy gaps between the ground and excited states. The energy gap correlation function can then be used to compute the spectral density, which shows how vibrational modes couple to the optical excitation and the resulting optical spectra, such as the resonance Raman spectrum and the absorption spectrum. This trajectory-based approach is able to capture both vibronic effects and specific environmental effects on the optical spectrum on an equal footing. Specific solvent interactions are known to cause spectral changes. The formation of hydrogen bonds (Hbonds) between chromophore and solvent can blue-shift or red-shift peaks in the optical spectrum. Recent studies have examined how H-bonding of water to flavin affects the vibrational frequencies and band intensities. One advantage of a trajectory-based approach is that changes in excitation energy and energy fluctuations can be connected to specific chromophoresolvent interactions and even solvent dynamics. For example, by careful analysis of the trajectory, for the anionic chromophore of photoactive yellow protein, deprotonated trans-thiophenyl-p-coumarate (pCT -), some of the authors showed that low-frequency motions of the chromophore with the surrounding solvent led to spectral broadening, and coupling of the high-frequency motion of the C=O stretch to solvent led to increased vibronic coupling. To fully realize the aforementioned advantages of the trajectory-based approach, we need an accurate description of both the chromophore and chromophore-solvent interactions, as well as a reliable approach for the computation of electronic excitation in the presence of solvent. For the latter, time-dependent density functional theory (TDDFT) has been the main workhorse with its modest computational cost and reasonable accuracy, and in our previous work, we have suggested guidelines for the use of TDDFT in excited state calculations of condensed phase systems. For the former issue, ab initio molecular dynamics (AIMD) with density functional theory (DFT) offers quantum-chemical accuracy at a high computational cost that limits the time scale for sampling configuration space, whereas more affordable force-field-based MD (FFMD), usually with empirical parameters, may suffer from poor descriptions of the chromophore or of specific chromophore-solvent interactions. Furthermore, it has been shown that a mismatch between the potential used for the MD simulation and the excitation energy calculations (e.g., using FFMD for the trajectory and TDDFT for excited state computation) introduces artifacts to the computed spectral density. It is thus desirable to either use the same quantum chemical protocol (i.e., functional, basis set, and solvent treatment) for MD simulation and excited state calculations (e.g. with AIMD and TDDFT employing the same functional), or to parametrize the force field (FF) using the same quantum chemical protocol as used for the excited state calculation to minimize the mismatch. Recent advances in ab initio based FF parametrization, such as the QUantum mechanical BEspoke Kit (QUBEkit), which derives system-specific FF parameters directly from quantum mechanical (QM) calculations, provide a promising way forward.