Many physical phenomena must be accounted for to accurately model solution-phase optical spectral line shapes, from the sampling of chromophore-solvent configurations to the electronic-vibrational transitions leading to vibronic fine structure. Here we thoroughly explore the role of nuclear quantum effects, direct and indirect solvent effects, and vibronic effects in the computation of the optical spectrum of the aqueously solvated anionic chromophores of green fluorescent protein and photoactive yellow protein. By analyzing the chromophore and solvent configurations, the distributions of vertical excitation energies, the absorption spectra computed within the ensemble approach, and the absorption spectra computed within the ensemble plus zero-temperature Franck-Condon approach, we show how solvent, nuclear quantum effects, and vibronic transitions alter the optical absorption spectra. We find that including nuclear quantum effects in the sampling of chromophore-solvent configurations using *ab initio* path integral molecular dynamics simulations leads to improved spectral shapes through three mechanisms. The three mechanisms that lead to line shape broadening and a better description of the high-energy tail are softening of heavy atom bonds in the chromophore that couple to the optically bright state, widening the distribution of vertical excitation energies from more diverse solvation environments, and redistributing spectral weight from the 0-0 vibronic transition to higher energy vibronic transitions when computing the Franck-Condon spectrum in a frozen solvent pocket. The absorption spectra computed using the combined ensemble plus zero-temperature Franck-Condon approach yield significant improvements in spectral shape and width compared to the spectra computed with the ensemble approach. Using the combined approach with configurations sampled from path integral molecular dynamics trajectories presents a significant step forward in accurately modeling the absorption spectra of aqueously solvated chromophores.

## I. INTRODUCTION

UV-Vis absorption spectra are sensitive to inter and intramolecular fluctuations, as well as to how the solvent environment couples with the chromophore’s ground and excited state energies. However, for systems that form strong hydrogen bonds, such as the green fluorescent protein (GFP) and photoactive yellow protein (PYP) chromophores, both of which are known to actively engage in hydrogen bonding in their native protein environments,^{1–4} continuum models of the solvation environment are insufficient. Modeling the absorption bands of chromophores in solution requires methods that go beyond only considering vertical excitation energies obtained from the minimum of the potential energy surface (PES).

One approach for modeling spectral shapes is to generate a temperature-dependent ensemble of configurations, e.g., from a molecular dynamics (MD) simulation. An excited state calculation is then performed on each configuration, yielding a distribution of vertical excitation energies. Each excitation can be convoluted with a Gaussian function to create an absorption spectrum that can then be compared to the experimental spectrum.^{5–9} This “ensemble approach” allows for the incorporation of a large number of explicit solvent molecules in the excited state electronic structure calculations. An explicit treatment of the electronic structure of large regions of solvent is known to be necessary for capturing solute-solvent polarization and charge-transfer interactions.^{7,9–17} Because this approach includes the fluctuating solvent environment, it provides temperature-dependent inhomogeneous broadening of the spectrum. However, for an ensemble approach calculation of the spectrum of the aqueous PYP chromophore, when the solute-solvent configurations are sampled using classical MD with a force field or classical *ab initio* MD (AIMD) with a density functional theory (DFT) treatment of the chromophore, the spectrum is too narrow^{10} and misses the high-energy shoulder observed experimentally.^{18,19}

However, nuclear quantum effects (NQEs) were not included in the configurational sampling in these simulations and these effects might account for some of the missing spectral width. Recent studies have included NQEs when generating the configurations used to compute absorption spectra within the ensemble approach. For example, NQEs were observed to give rise to a red shift and broadening in the spectra of systems ranging from gas-phase neat water clusters of varying size^{20–22} to biological systems like ellipticine,^{23} the DNA base methylguanine,^{24} and amyloid protein fibrils.^{25} A very recent study by Law and Hassanali^{26} showed that the amount of spectral broadening due to NQEs was similar for chromophores with different hydrogen bond donating and accepting abilities.

To the best of our knowledge, previous studies of the role of NQEs on optical absorption spectra have used the ensemble approach to obtain a ground state distribution of vertical excitation energies from solute-solvent configurations. This approach neglects vibronic transitions and thus will not capture any vibronic fine structure in the optical spectrum. Our recently developed ensemble-zero temperature Franck Condon (E-ZTFC) approach for modeling spectral shapes combines the sampling of an ensemble of solute-explicit solvent configurations with a zero-temperature Franck-Condon shape function.^{27} This approach accounts for both the temperature-dependent inhomogeneous broadening due to solute-solvent interactions and the vibronic transitions that broaden the spectrum and lead to spectral asymmetry.

Here, we analyze how NQEs, both in the sampling of the ground state nuclear configurations and in treating the vibrational transitions accompanying the electronic transitions, change the optical absorption spectral line shape of aqueously solvated GFP and PYP anionic chromophores (Fig. 1). To efficiently sample nuclear configurations including NQEs, we use our ring polymer contraction (RPC) approach^{28,29} to perform *ab initio* path integral molecular dynamics (AI-PIMD) simulations^{30} of the chromophores both in water and vacuum. By comparing the chromophore structures and the distributions of excitation energies of the solvated and vacuum trajectories, we then assess both the direct and indirect influence of the solvent environment in the form of solvent polarization effects and solvent-induced conformational changes of the chromophore. We then apply our recently introduced E-ZTFC approach for capturing absorption lineshapes^{27} to analyze the importance of vibronic transitions when modeling the spectra of aqueously solvated chromophores. Overall, we find that NQEs lead to spectral broadening in three ways: by softening the bonds between the heavy atoms of the chromophores, by giving rise to a wider range of vertical excitation energies due to direct solvent interactions, and by shifting the intensity of the vibronic transitions accompanying each vertical excitation.

## II. METHODS AND COMPUTATIONAL DETAILS

### A. Calculation of the absorption spectrum within the combined ensemble plus zero-temperature Franck-Condon approach

Within the E-ZTFC approach, solute-solvent configurations are sampled just as in the ensemble approach, i.e., by performing simulations at the desired temperature using a fully flexible solute and the surrounding solvent. Here we sample configurations via AI-PIMD, which includes NQEs such as zero-point energy and also allows for bond making and breaking as dictated by the electronic structure method. Such an approach should capture proton sharing and transfer between the solute and solvent.

Performing excited state calculations for each configuration and convoluting each excitation with a Gaussian function yields an optical spectrum computed with the traditional ensemble approach. This approach neglects vibronic transitions. Here we also compute spectra using our recently introduced E-ZTFC approach. In this approach, instead of Gaussian broadening, we account for vibronic fine structure by convoluting a zero-temperature Franck-Condon (ZTFC) shape function with the excitation energies computed from an ensemble of solute-solvent configurations obtained from an MD trajectory at the desired temperature. The Franck-Condon spectrum at zero-temperature is given by

Here, *μ*_{if} is the electronic transition dipole moment of the transition going from initial electronic state *i* to final electronic state *f*, and $\varphi vi$ and $\varphi vf$ denote nuclear wave functions of the Born-Oppenheimer PES of electronic states *i* and *f*, respectively. $Evf$ denotes the total energy of the system in electronic state *f* and vibrational mode $vf$. Since the vibronic spectrum is computed in the zero-temperature limit, there is a sum over vibrational modes of the final state only because all transitions arise from the ground vibrational state. The symbol $N$ denotes a Gaussian function with mean $Evf\u2212Evi0$ and standard deviation *σ*. In the E-ZTFC approach, all the temperature effects are accounted for via an integral over solute-solvent configurations. To compute the absorption spectrum, we first rewrite the average excitation energy $Evf\u2212Evi0$ for a given Franck-Condon transition as^{27}

where $E0i\u2192vfvib$ denotes the purely vibrational contribution to the vibronic energy in transitioning from the ground state vibrational mode of the ground state PES to the vibrational mode $vf$ of the excited state PES, Δ_{if} is the adiabatic excitation energy, and $Efluctav(T)$ is the average energy due to temperature fluctuations sampled through the conformational integral. *ω*_{if} denotes the vertical excitation energy of the electronic transition going from initial state *i* to final state *f*, **R** is a collective variable for all solute and solvent nuclear degrees of freedom, and *ρ*_{GS}(**R**, *T*) denotes a probability distribution function for a given solute-solvent conformation **R**, with the electronic system in its ground state, at temperature *T*. $R0GS$ is the optimized ground state structure on the ground state PES.

Temperature-induced fluctuations are included in the computed absorption spectrum by sampling the conformational integral using MD performed at the desired temperature. Extracting a set of *N*_{frames} independent solute-solvent conformations {**R**_{j}} from an MD trajectory, the full absorption spectrum in the E-ZTFC approach can then be written as^{27}

where *f*_{if} is the electronic oscillator strength of transition *i* → *f* and *σ* is a small numerical convergence parameter that accounts for finite sampling. The final expression reduces to a convolution of a single ZTFC spectrum with a set of vertical excitation energies computed for an ensemble of conformations sampled with MD, where the solute and solvent degrees of freedom are fully coupled. The E-ZTFC approach thus is not much more computationally expensive than the ensemble approach alone, requiring only the additional computational cost of the Franck-Condon spectrum (or spectra, if averaging over multiple Franck-Condon spectra as is done here) used to construct the vibronic shape function.

The E-ZTFC approach as outlined above relies on the Franck-Condon approximation being valid for computing the zero-temperature vibronic shape function and is thus appropriate for bright transitions. In principle, the approach can be expanded to weakly allowed and dipole-forbidden states by including Herzberg-Teller effects^{31} in Eq. (1). However, for the purpose of this work, the focus is on the bright S_{1} transition of the GFP and PYP chromophore anion, where the expression in Eq. (1) is sufficient to describe the vibronic fine structure of the transition.

Computing absorption spectra within the E-ZTFC approach therefore requires three steps. The first is the generation of the conformations of the chromophore in the solvent for a given temperature. In this work, we obtain the configurations using both AIMD and AI-PIMD in order to assess the influence of ground state NQEs on the computed spectra. The second is the calculation of vertical excitation energies and oscillator strengths for the solute-solvent conformations. The third is the computation of a vibronic shape function. In Sec. II B, we first examine issues related to double counting of the nuclear degrees of freedom of the chromophore in the E-ZTFC approach and then provide the computational details for each of the three steps in Secs. II C, II D, and II E.

### B. Double-counting of nuclear degrees of freedom of the chromophore in the E-ZTFC approach

The E-ZTFC approach includes chromophore nuclear degrees of freedom both from sampling the chromophore motion in the ensemble of vertical excitation energies and from the ground state vibrational wave function used to construct the zero temperature Franck-Condon shape function. To characterize how this double-counting of the chromophore nuclear motion affects the computed spectra, we introduced a simple model system in Ref. 27, consisting of an ensemble of identical displaced harmonic oscillators coupled to a uniform classical solvent environment described by the solvent reorganization energy. We found that when the ensemble spectrum was sampled from a Boltzmann distribution, the double-counting in the E-ZTFC approach led to systematic overestimation of the spectral width compared to the exact solution of the model system but that this overestimation was small in the limit of large solvent reorganization energy.

In the present work, the ensemble of vertical excitation energies is sampled from both AIMD and AI-PIMD trajectories. This AI-PIMD distribution includes nuclear quantum effects, including zero-point energy, allowing the system to sample more anharmonic regions of the potential energy surface. To test the validity of combining the E-ZTFC approach with the AI-PIMD ensemble that accounts for nuclear quantum effects on the ground state potential energy surface, here we again explore the effects of double-counting the nuclear motion of the chromophore. We extend our simple model system to include a quantum distribution of the nuclear degrees of freedom. Furthermore, to assess the behavior of the E-ZTFC approach for anharmonic degrees of freedom, we also extend the model system to that of an ensemble of identical displaced Morse oscillators. In Sec. I of the supplementary material, we provide the mathematical details of this model, results for a variety of solute-solvent interaction strengths, and additionally compare computed spectral results for the GFP chromophore with and without double-counting for a subset of configurations. Furthermore, in Sec. I B of the supplementary material, we show in detail how the E-ZTFC approach is an approximation to the computationally much more demanding, double-counting free approach that can be obtained when applying a separation of time scales to solute and solvent motion. Time scale separation arguments have previously been applied to account for the slow, anharmonic degrees of freedom of flexible chromophores through ensemble sampling.^{32,33} Here, we extend this argument to separate the time-evolution of the entire solvent environment from the vibrations of the chromophore such that the instantaneous excitation of the chromophore can be computed in the static field provided by a given solvent configuration.

Here, we consider a model consisting of an ensemble of oscillators with identical shape and displacement, for both harmonic and anharmonic potential energy surfaces, coupled to a uniform solvent environment. Two different forms of solvent broadening are considered in the model system. The first is broadening described purely by solvent electrostatic fluctuations through the solvent reorganization energy *λ*; this is the dominant broadening expected for solvents interacting weakly with the chromophore where the nuclear motion is fully independent from solvent degrees of freedom. For the second model of the solvent environment, we add an additional term *δ* to the solvent broadening meant to describe the direct interaction between the solute and the solvent. A model system defined in this way can be related to the more rigorous double-counting free approach that is considered in detail in Secs. I B and I C of the supplementary material for the GFP chromophore. The spectra for the model systems under the two different solvent models can be obtained exactly for both the harmonic and the Morse oscillator and can be compared to the corresponding E-ZTFC approach.

The simulated spectra for both harmonic and anharmonic model systems are shown in Fig. 2, which shows the exact solution for these model systems and the E-ZTFC spectra based on a quantum and a classical ensemble distribution. For the harmonic oscillator in Fig. 2(a), the finite temperature Franck-Condon (FTFC) result is the exact result, whereas for the Morse oscillators in Figs. 2(b) and 2(c), the FTFC spectrum is determined within the harmonic approximation to the ground and excited state potential energy surfaces, resulting in spectra that are too narrow compared to the exact solutions. The spectra in Figs. 2(a) and 2(b) include solvent broadening due to the solvent reorganization energy *λ*. In Fig. 2(c), for the exact results and the E-ZTFC spectra, we also include the effect of direct solute-solvent interaction through the additional broadening term (see Sec. I A of the supplementary material for further details and Sec. I C of the supplementary material where we provide a physical justification for the values of the solvent broadening terms based on the results for the GFP chromophore). In Fig. 2(c), we also show the FTFC spectrum where the solvent broadening originates only from the electrostatic contribution *λ*, ignoring any direct solute-solvent interactions. The results show that in Figs. 2(a) and 2(b), where the motion of the chromophore is considered to be fully independent from the fluctuations of the solvent environment, both E-ZTFC approaches systematically overestimate the spectral width, with the overestimation being more pronounced in the quantum ensemble sampling due to double-counting of zero point motion. For the anharmonic Morse oscillator, Fig. 2(b), the E-ZTFC approach based on the classical ensemble partially captures the high energy tail missing from the Franck-Condon spectrum based on the harmonic approximation. However, once we account for the effect of direct solute-solvent interactions in our model system, Fig. 2(c), we find that the E-ZTFC approach based on the quantum ensemble only slightly overestimates the full width of the spectrum and the classical E-ZTFC approach underestimates the spectral width. Thus, in the anharmonic system, the quantum ensemble broadening allows the E-ZTFC approach to approximately capture the correct anharmonic high energy behavior.

The systems studied in this work are anionic chromophores in water with strong H-bonding interactions and thus can be expected to have strong solute-solvent interactions. In Sec. I C of the supplementary material, we show that the model system including direct solute-solvent effects more accurately reflects the results for the chromophores studied here compared to the model system with only electrostatic broadening. Thus, even though the AI-PIMD sampling leads to increased double-counting of the nuclear degrees of freedom of the chromophore in high frequency modes due to the inclusion of the zero-point motion both in the ensemble sampling and the harmonic vibronic shape function, the AI-PIMD ensemble sampling partially captures the effect of anharmonicities not included in the harmonic ZTFC spectrum and only partially included from sampling from AIMD snapshots. Also, the AI-PIMD trajectories more accurately take into account specific chromophore-solvent interactions and sampling of anharmonic regions of the potential energy surface, both of which contribute to the spectra in the E-ZTFC approach.

In general, we expect the E-ZTFC approach based on AI-PIMD sampling of the nuclear degrees of freedom to yield good results for semi-flexible systems coupling strongly to their solvent environment. In more weakly coupled systems, where the motion of the chromophore is not as strongly influenced by the presence of the solvent, the E-ZTFC approach presented here might yield a more considerable overestimation of spectral broadening. In these cases, a rigorously double-counting free approach (see Sec. I B of the supplementary material), albeit computationally considerably more expensive, may be necessary to yield more reliable results.

### C. *Ab initio* and path integral molecular dynamics simulations

We performed classical and path integral AIMD simulations of the two anionic chromophores, both in vacuum and in water in the NVT ensemble at *T* = 300 K. For the aqueous simulations, periodic boundary conditions were used and the electronic structure of the system was described using DFT. Simulations were performed using the i-PI program^{34} and employed a multiple time scale (MTS) integrator of the r-RESPA form.^{35}

Classical AIMD MTS simulations employed a 2.0 fs outer time step and a 0.5 fs inner time step. AI-PIMD simulations employed ring polymer contraction (RPC)^{28–30} with centroid contraction, i.e., contraction to *P*′ = 1 replicas of the system and an MTS propagator with an outer time step of 2.0 fs and an inner time step of 0.25 fs. PIMD simulations^{36–38} were performed by thermostatting the non-centroid normal modes using target Langevin thermostats within the path integral Langevin equation (PILE) approach.^{38}

In the RPC and MTS simulations, the full forces were evaluated using the CP2K program^{39,40} and the revPBE density functional,^{41,42} with D3 dispersion corrections^{43} added. Atomic cores were represented using the dual-space Goedecker-Teter-Hutter pseudopotentials.^{44} The GPW method^{45} used Kohn-Sham orbitals expanded in the TZV2P basis set, and an auxiliary plane-wave basis with a cutoff of 400 Ry was used to represent the density. The criterion for convergence of the self-consistent field was an electronic gradient below a tolerance of *ϵ*_{SCF} = 5 × 10^{−7} using the orbital transformation method,^{46} with the initial guess provided by the always-stable predictor-corrector extrapolation method^{47,48} at each MD step. The RPC and MTS reference forces were evaluated at the SCC-DFTB3^{49} level of theory using the DFTB+ program.^{50} The 3ob parameter set^{51} was used and dispersion was included via a Lennard-Jones potential^{52} with parameters taken from the Universal Force Field.^{53}

Initial configurations were obtained from classical molecular mechanics (MM) simulations performed with the Sander module of the AmberTools16 software package.^{54} Force field parameters for the chromophores were obtained from the general Amber force field,^{55,56} using the force field for the PYP chromophore described in previous work^{10} and the reparameterization for the GFP chromophore described in Ref. 27. The q-SPC/Fw water model^{57} was used. The PYP chromophore was solvated in a rectangular box by 166 water molecules and the GFP chromophore was solvated in a rectangular box by 167 water molecules. Classical MM simulations of 20 ns in length were performed in the NPT ensemble at *T* = 300 K in order to compute the average densities. In these MM simulations, chromophores were restrained to lie along the longest axis of the solvent box via harmonic biasing potentials using the PLUMED software package.^{58} Additional classical MM simulations were then performed in the NVT ensemble at *T* = 300 K at the average densities obtained from the NPT simulations in order to obtain 10 uncorrelated initial configurations for each chromophore that were spaced by 1 ns in time.

Using DFT-based classical AIMD, each configuration was equilibrated for 3 ps using a local Langevin thermostat with a time constant of 25 fs. Production runs were then performed using a local Langevin thermostat with a time constant of 2 ps. For each chromophore, 10 classical AIMD and 10 AI-PIMD simulations of 15 ps in length were performed, with the AI-PIMD simulations launched from the final configurations of the classical AIMD simulations. In aggregate, 300 ps of classical AIMD and 300 ps of AI-PIMD were performed for the two aqueously solvated chromophores.

In order to obtain the final configurations used to compute the absorption spectra, we first extracted 100 evenly spaced configurations from each of the 15 ps trajectory segments, for a total of 1000 snapshots for AIMD and 1000 snapshots for AI-PIMD for each chromophore. The chromophore was centered in the periodic box and any solvent water molecules that were broken across the edges of the periodic box were reconstructed. Because we wanted to include large amounts of explicit solvent in the calculations of the vertical excitation energies, we added more water molecules beyond the ∼166 originally solvating the chromophores during the MD. In order to increase the amount of solvent around the chromophore, each DFT solute-solvent configuration was centered in a larger cubic box and solvated with additional ∼6000-7000 q-SPC/Fw water molecules such that the edges of the original DFT box were at least 20 Å from the edges of the larger cubic box. One of the box vectors of this larger box was then elongated by a factor of 3 to provide a vacuum region of the cell into which the water could expand if needed (refer to Sec. II of the supplementary material). This system containing the additional solvent was then equilibrated in the NVT ensemble at T = 300 K for 1 ns, with the original centered DFT solute-solvent configuration held frozen. We utilized the OpenMM^{59–62} software package to perform the relaxation of this additional solvent. For this simulation, we employed a Langevin integrator with a friction coefficient of 1 ps^{−1}. The centered frozen DFT configuration remained at the center of the open slab over the course of the equilibration. The final configurations obtained from the end of these 1000 relaxation runs were then used to calculate the absorption spectra. Validation of this procedure regarding the convergence of the absorption spectrum with regard to the DFT box size and amount of MM solvent is discussed in the supplementary material. Furthermore, Sec. IV of the supplementary material contains a detailed analysis of the convergence of vertical absorption spectra with respect to the 10 individual uncorrelated AIMD and AI-PIMD trajectories, as well as a comparison of the AIMD ensemble spectrum to the ensemble spectrum obtained from a classical MM trajectory of 8 ns^{27} for the GFP chromophore anion, demonstrating the robustness of the sampling protocol chosen in this work.

To assess the indirect effects of the solvent on the configurations sampled by the chromophore, gas-phase simulations of the chromophores were also performed. These calculations used the i-PI program^{34} and employed a multiple time scale (MTS) integrator of the r-RESPA form.^{35} For these simulations, the DFT energies and forces were computed using the TeraChem software package.^{63} We employed the revPBE density functional,^{41,42} with D3 dispersion corrections^{43} added. Atom-centered orbitals were expanded in the 6-311++G(d,p) basis set, and a wave function convergence threshold of *ϵ* = 1 × 10^{−6} was used. Classical simulations employed an outer time step of 2.0 fs, an inner time step of 0.5 fs, and a local Langevin thermostat with a time constant of 0.5 ps. AI-PIMD simulations employed an outer time step of 2.0 fs and an inner time step of 0.25 fs, and non-centroid normal modes were thermostatted using target Langevin thermostats within the path integral Langevin equation (PILE) approach.^{38} We accumulated 100 ps of both AIMD and AI-PIMD for each chromophore, discarded 5 ps of each trajectory for equilibration, and sampled frames using a stride of 150 fs (75 frames) for the analysis in order to maintain consistency with the analysis of the condensed phase systems.

Experimentally, the PYP anion chromophore is predominantly protonated at pH = 7, whereas for the GFP chromophore, both the anionic and the neutral form have appreciable populations in equilibrium at pH = 7. Because proton transfer from the water to the chromophore anions did not occur in our simulations, we therefore compare our computed spectra to experimental spectra obtained in an aqueous solution with a borate buffer at pH = 10.2 for the PYP chromophore^{19} and an aqueous solution with NaOH at pH = 13 for the GFP chromophore.^{64} These experimental conditions guarantee that only the anion is present in solution.

### D. Vertical excitation energies

The direct solute-solvent interactions were accounted for in the computation of the excitation energies by explicitly including a significant part of the solvent environment in the QM region of the electronic structure calculation, with the remaining solvent included as MM fixed point charges. The size of the QM region was controlled by the cutoff radius *R*_{cut} = 8 Å, where all solvent molecules with a center of mass within *R*_{cut} of any atom of the chromophore were treated with QM. All solvent molecules with centers of mass outside of the cutoff radius were included in the excited state calculation as electrostatic point charges.

The lowest three time-dependent density functional theory (TDDFT) singlet states were computed for all snapshots; three states was found to be sufficient to cover the entire range of the visible spectrum. All vertical excitation energy calculations were carried out using linear response TDDFT as implemented in the TeraChem software package.^{65} The Tamm-Dancoff approximation^{66} was used throughout and all calculations were performed with the range-separated hybrid density functional CAM-B3LYP^{67} and the 6-31+G(d) basis set. The default range-separation parameter 0.33 was used for all calculations. As previously noted,^{68} the basis set is too small to fully converge the absolute values of the excitation energies in these systems, which is especially apparent for the GFP anion. However, the present work focuses on the shape of absorption spectra rather than the absolute position of absorption maxima; the spectral shape has been shown to be relatively robust with respect to changes in the basis set size.^{27}

### E. Vibronic shape function

The explicit QM solvent was frozen while computing the Franck-Condon spectra for a number of solute-solvent conformations corresponding to an approximation of a separation of time scales between solute vibrations and solvent relaxation effects. Denoting the frozen solvent coordinates of a given snapshot **R**_{j} as $Rj{m}$, an average normalized vibronic shape function $\sigma ZTFCav$ can be defined via

where *c*_{norm} is a normalization constant chosen such that $\sigma ZTFCav$ integrates to 1. In Eq. (4), the individual Franck-Condon spectra *σ*_{ZTFC} are now computed in different effective frozen solvent environments. The specific frozen solvent environment of a given snapshot changes the equilibrium geometry of the chromophore, as well as the chromophore vibrational wave functions and energies. This influence of a specific solvent configuration on the computed shape function is denoted by the parametric dependence on $Rj{m}$ that has been introduced in the expression for *σ*_{ZTFC} in Eq. (4). $E0\u22120Rj{m}$ is the total energy of the 0-0 transition for the Franck-Condon spectrum calculated in frozen solvent environment $Rj{m}$ and $E0\u22120av$ is the average 0-0 transition energy for the *N*_{shape} frozen solvent snapshots. The spectra are averaged by making their *E*_{0−0} transition energies coincide, corresponding to the assumption that variations in the *E*_{0−0} transition energy for different solvent environments are accounted for in the variations of the TDDFT vertical transition energies. The resulting average normalized vibronic shape function is then used in Eq. (3) such that the final expression for the E-ZTFC computed absorption spectrum reduces to

where $\omega ifav=1Nframes\u2211jNframes\omega if(Rj)$. Since $\sigma ZTFCav$ is normalized to 1, all contributions to the oscillator strength of the convoluted spectrum arise from the oscillator strengths $fifRj$ of the individual ensemble snapshots. Thus the E-ZTFC approach takes into account variations in the oscillator strength due to different solute-solvent configurations such that the Condon approximation of a constant electronic transition dipole moment is only applied to the chromophore degrees of freedom when calculating the individual zero temperature Franck-Condon shape functions.

In this work, the vibronic shape function was computed according to Eq. (4) as an average of five explicit solute-solvent conformations chosen from the uncorrelated MD snapshots. Using only five frozen solvent environments will not necessarily generate a shape function that represents the true average vibronic fine structure of the system, especially if different solvent environments can produce significant fluctuations in the vibronic fine structure, as is likely the case in systems with strong solute-solvent interactions.^{27} However, in these types of systems, where the vibronic spectrum of the solute embedded in the solvent pocket is strongly dependent on the solvent configuration, the use of a single average vibronic shape function for all configurations sampled in the E-ZTFC approach is likely the most dominant error introduced in the calculations, rather than errors in the computed average shape function itself due to a lack of sampling.

The amount of QM solvent included in the calculation of the vibronic shape function was defined through the cutoff radius *R*_{cut} in the same way as for the computation of vertical excitation energies. However, because of the high computational cost associated with optimizing the ground and excited state geometry of the chromophore, a cutoff radius of *R*_{cut} = 3 Å, which is significantly smaller than the cutoff radius for the vertical excitation energies, was used for the frozen solvent region during the geometry optimizations. This explicit solvent representation included the solvent molecules that interact most strongly with the solute, such as hydrogen-bonded solvent molecules, but was not large enough to provide a full solvation shell. The QM solute and solvent region was surrounded by a polarizable continuum model.^{69}

Geometry optimizations and frequencies, as well as zero temperature Franck-Condon spectra, were computed using the Gaussian development version^{70} at the CAM-B3LYP/6-31+G(d) level of theory, consistent with the computation of vertical absorption energies. Only the vibrational frequencies of the chromophore were considered in the Franck-Condon spectra. Since the chromophore degrees of freedom were optimized in the frozen solvent pocket, this allowed the rigorous definition of the harmonic frequencies of the chromophore in terms of the Hessian around the ground and excited state optimized geometries. As in the Terachem calculations of vertical excitation energies, the range-separation parameter was 0.33 for the CAM-B3LYP functional. During the geometry optimization and frequency calculations, the QM solvent atoms were kept frozen and dispersion interactions between the molecules were accounted for using Grimme’s D3 empirical dispersion correction.^{43} All Franck-Condon spectra were broadened by a Gaussian function with standard deviation *σ* = 0.0105 eV. This broadening parameter, together with a sampling of vertical excitation energies for 1000 snapshots, was found to be sufficient to yield smooth absorption spectra, while preserving all of the dominant features in the vibronic fine structure of the Franck-Condon spectra. For both chromophores, the TDDFT S_{1} excited state was the only state with significant oscillator strength and was thus the only state convoluted with the vibronic shape function. A Gaussian function with *σ* = 0.0105 eV was used for all other excitations.

The quality of the computed Franck-Condon spectra was assessed by analyzing the sum of the individual Franck-Condon intensities in comparison to the total expected intensity following analytic sum rules. Average spectrum progressions of 96.5% and 91.6% were found for the GFP chromophore using AIMD and AI-PIMD snapshots, respectively. For the PYP chromophore, the average spectrum progressions were 93.8% and 96.8% for the AIMD and AI-PIMD snapshots, respectively.

## III. RESULTS AND DISCUSSION

We first present analysis of the vertical excitation energies obtained from computations on snapshots from the AIMD and AI-PIMD trajectories. Second, we focus on understanding the solvent-induced changes in the excitation energies. To simplify our analysis, in Secs. III A and III B, we compute spectra within the ensemble approach, not including any vibronic effects due to the simultaneous excitation of electronic and vibrational states. Third, in Sec. III C, we compare the zero-temperature Franck-Condon vibronic spectra computed from configurations extracted from AIMD and AI-PIMD dynamics. Finally, in Sec. III D, we construct the full vibronically broadened absorption spectra in the combined E-ZTFC approach for the AIMD and AI-PIMD data sets for both chromophores and compare the computed spectra to experimental results.

### A. The role of NQEs on the vertical excitation energies

We calculated the vertical excitation energies for configurations extracted from our AIMD and AI-PIMD simulations of the PYP and GFP chromophores in vacuum and water. Figure 3 shows the vertical excitation energies and the corresponding absorption spectra in water obtained using the ensemble approach, wherein the excitation energies are broadened by a Gaussian with a standard deviation of *σ* = 0.0105 eV and the height of the Gaussian is determined by the value of the oscillator strength of the transition. For the aqueously solvated PYP and GFP chromophores, NQEs lead to a red-shifting of the average excitation energy by 0.073 eV and 0.050 eV, respectively. This red shift in solution due to NQEs is similar to that obtained in vacuum (0.058 eV and 0.068 eV for the PYP and GFP chromophores, respectively). The shift in the average excitation energy is also accompanied by significant broadening and increased asymmetry of the spectrum, as well as a more prominent high-energy shoulder (see Sec. VI of the supplementary material for a detailed analysis of non-Gaussian features in the distribution of vertical excitation energies). To quantify the broadening due to NQEs and solvation, Table I shows the changes in the standard deviation of the excitation energy distributions relative to the vacuum value. The increase in the standard deviations of the distributions by ∼0.03 eV upon quantizing the nuclei is consistent for both PYP and GFP chromophores across all environments. In the case of the GFP chromophore, the increase in the width arising from the inclusion of NQEs is as large as that obtained from solvating the molecule and hence is not a small perturbation, whereas for PYP the solvation effect is ∼2 fold larger than that due to NQEs. The larger degree of inhomogeneous broadening for the PYP chromophore compared to the GFP chromophore is likely due to a larger rearrangement of charge upon electronic excitation, which is evidenced by a larger average magnitude for the transition dipole moment and a larger difference in the ground and excited state dipole moments.

. | PYP . | GFP . | ||
---|---|---|---|---|

. | AIMD . | PI-AIMD . | AIMD . | AI-PIMD . |

Vacuum | 0.000 | 0.031 | 0.000 | 0.035 |

Solvated | 0.071 | 0.106 | 0.039 | 0.096 |

Stripped | 0.016 | 0.046 | 0.007 | 0.052 |

Solvent | 0.042 | 0.062 | 0.011 | 0.051 |

. | PYP . | GFP . | ||
---|---|---|---|---|

. | AIMD . | PI-AIMD . | AIMD . | AI-PIMD . |

Vacuum | 0.000 | 0.031 | 0.000 | 0.035 |

Solvated | 0.071 | 0.106 | 0.039 | 0.096 |

Stripped | 0.016 | 0.046 | 0.007 | 0.052 |

Solvent | 0.042 | 0.062 | 0.011 | 0.051 |

To elucidate the structural origins of these changes, we therefore assessed which degrees of freedom give rise to the largest changes in the vertical excitation energies. For both the PYP and GFP chromophores, we found that the structural parameter that most strongly correlates with the excitation energies is the length of the C—C bond connecting the phenolate ring to the conjugated backbone, labeled C1—C2 in Fig. 1. Furthermore, for the GFP chromophore, there is experimental evidence that the C—N double bond in the imidazole, labeled C3—N1 in Fig. 1(b), couples strongly to the bright state^{71} and it is therefore included in our analysis. The comparable increase in the standard deviation of excitation energies in both the vacuum and solvated systems upon including NQEs observed in Table I is consistent with the structural parameters of the chromophore, rather than the solvent, being strongly coupled to the bright transition. To simplify our analysis of the correlation of these bond lengths with the value of the excitation energy, we removed the inhomogeneous solvent environment, resulting in “stripped” chromophore configurations, and then recomputed the excitation energies.

Figures 4(a) and 4(b) show the correlation between the vertical excitation energy of the bright S_{1} state and the C—C bond length for the PYP chromophore and the combined C—C and C—N bond lengths for the GFP chromophore, respectively. In both cases, increasing the bond lengths leads to an decrease in the excitation energies, and the slopes of the best fit lines are almost identical for the AIMD and AI-PIMD simulations. From these contour plots, it is immediately clear that the AI-PIMD simulations sample a significantly broader range of bond lengths. Although these are bonds between relatively heavy atoms, the zero-point energy added when NQEs are included is significant, ∼1.7 kcal/mol for each C—C stretch and ∼2.6 kcal/mol for the C—N double bond. The zero-point energy included in the AI-PIMD simulations therefore leads to an increase in the breadth of the bond length distributions of ∼50% for both the C—C bond connecting the phenolate to the conjugated backbone and the C—N double bond in the imidazole in the GFP chromophore (see Sec. VII of the supplementary material). The increased range of bond lengths of the solute in the AI-PIMD simulations is thus a major source of the spectral broadening when including NQEs.

It is worth noting that we also considered other chromophore structural parameters that might be expected to correlate with the excitation energies. Examples include the dihedral angle connecting the phenolate to the conjugated backbone, because strong twists have the potential to break conjugation and thus alter the character of the S_{1} state, as well as the C—O bond length of the phenolate oxygen. However, although we found that both NQEs and solvation alter these degrees of freedom (see Sec. VII of the supplementary material), no correlation with respect to excitation energies was obtained.

### B. The role of the solvent environment on the vertical excitation energies

We now consider how the solvent environment changes the vertical excitation energies. For both chromophores, the addition of solvent blue-shifts the mean excitation energy: for PYP by 0.083 eV and 0.068 eV in AIMD simulations and AI-PIMD simulations, respectively, and for GFP by 0.031 eV and 0.049 eV. As shown in Table I, the breadth of the excitation energy distributions also increases markedly upon solvation. To assess the origins of these changes, we consider the excitation energy distributions obtained from configurations of the solvated chromophore with the solvent “stripped” away, leaving only the chromophore. Because the internal degrees of freedom of the chromophores are unchanged between the snapshots, the differences in vertical excitation energies must arise from the interactions of the solvent with the chromophore. Hence these “stripped” configurations allow us to assess whether changes in the computed excitation energies arise indirectly from changes in the structure of the chromophore or directly from chromophore-solvent interactions. The computed ensemble spectra for the stripped configurations and for the fully solvated configurations are shown in Fig. 5, and the computed standard deviations of the distribution of vertical excitation energies are given in Table I.

The most notable difference between the solvated and stripped spectra in Fig. 5 is the increased width for the solvated spectra. The high-energy shoulder in the solvated spectra is missing in the stripped spectra, suggesting that this spectral feature originates from direct chromophore-solvent interactions. To quantify the indirect and direct effects of solvent on the broadening, we analyze the standard deviations of the excitation energy distributions in more detail.

Table I shows that the stripped configurations have a small increase in the standard deviation of excitation energies compared to the vacuum ones (see Sec. VIII of the supplementary material for a more detailed comparison of the ensemble spectra computed for the vacuum and the stripped snapshots). The similar standard deviation for the stripped and vacuum excitation energies indicates that the solvent’s ability to stabilize configurations of the chromophore that give rise to a wider range of excitation energies is relatively mild. For both chromophores, there is an increase in the number of conformations where the phenolate is twisted with respect to the conjugated backbone for the stripped configurations over the vacuum ones (see Sec. V of the supplementary material), but the observed changes do not significantly alter the character of the S_{1} state.

The majority of the broadening of the distribution of vertical excitation energies upon solvation therefore arises from the direct effect of the solvent on the chromophores’ electronic structure. By comparing the spread of vertical excitation energies for the solvated snapshots and for the same snapshots with the solvent stripped away, we can estimate the broadening due to the direct solvent environment. The standard deviation of the excitation energies due to the solvent environment can be quantified by approximating the solvated distribution as a convolution of the stripped spectrum with a Gaussian to account for the solvent broadening. Using the well-known result for the convolution of two Gaussians, under the above approximations, we can write the standard deviation *σ*_{solvent} of the solvent broadening as

where *σ*_{solvated} and *σ*_{stripped} correspond to the standard deviations of the distribution of vertical excitation energies for the solvated snapshots and the same snapshots where the solvent has been stripped away. Using the values of *σ*_{solvated} and *σ*_{stripped} from Table I, the estimated effective solvent broadening for both chromophores shows a significant increase in the effective solvent broadening when accounting for NQEs, with AI-PIMD *σ*_{solvent} increasing by 15% compared to the AIMD results for the PYP chromophore and by 34% for the GFP chromophore. This increase suggests that the broadening of the spread of vertical excitation energies when accounting for NQEs is not only due to the internal motion of the solute, specifically the softening of C—C and C—N bonds as discussed in Sec. III A, but also due to the solvent. Thus, although the net blue shift of the spectrum caused by solvent effects is similar for the AIMD and AI-PIMD data sets, the solvent-induced broadening is stronger in the AI-PIMD configurations.

To assess the contribution of hydrogen bonding to the solvent-induced spectral broadening, we looked at the hydrogen bonds from the solvent to the phenolate oxygen. As shown in Sec. VIII of the supplementary material, the NQEs increase the proton sharing between the phenolate and solvent, but the proton-sharing coordinate does not correlate strongly with the excitation energy. This lack of correlation suggests that although NQEs yield a significant increase in solvent-induced broadening in both systems, this increase is most likely a collective effect due to the structure of the first solvation shell and cannot be attributed to a single hydrogen bonding site.

### C. The role of vibronic transitions

In addition to an ensemble of vertical excitation energies, the calculation of a zero-temperature Franck-Condon (ZTFC) shape function is necessary for computing spectra within the E-ZTFC approach. In this work, the ZTFC shape function is computed as an average over five snapshots taken from uncorrelated MD trajectories. The average vibronic shape functions are reported in Fig. 6(a) for the PYP chromophore anion and Fig. 6(b) for the GFP chromophore anion. The individual ZTFC spectrum for each snapshot is reported in Sec. IX of the supplementary material.

The spectra reported in the supplementary material demonstrate that there is significant variation between individual ZTFC spectra for different frozen solvent environments, indicating that the specific solvent environment has a large influence on vibronic transitions in the solutes. Using a single vibronic shape function computed as an average from the Franck-Condon spectra in different environments is therefore certainly an approximation in the E-ZTFC approach but one that leads to significant computational savings over computing a Franck-Condon spectrum for each snapshot in the ensemble. For the GFP chromophore anion, the average ZTFC spectra for the AIMD and AI-PIMD trajectories agree remarkably well. For the PYP chromophore anion, however, the AI-PIMD vibronic shape function distributes more spectral weight from the first 0-0 vibronic peak to higher energy transitions compared to the AIMD vibronic shape function. The structural changes between the ground and S_{1} state optimized geometries are related to the strength of the hydrogen bonding of water molecules to the phenolate oxygen. In configurations with strong hydrogen bonds in the ground state, these hydrogen bonds become less favorable in the excited state, causing a larger reorientation of the chromophore (see Sec. IX of the supplementary material).

The shift in spectral weight from the 0-0 transition to higher energy vibronic states for the AI-PIMD shape function of PYP broadens the absorption spectrum computed within the E-ZTFC approach (see Sec. III D). For the spectrum computed from AI-PIMD snapshots, there will therefore be an increase in width compared to the spectrum computed with AIMD snapshots for three reasons: (1) from the increased spread in vertical excitation energies due to the spread in heavy atom bond lengths of the chromophore, (2) from the increased spread in vertical excitation energies due to the collective effect of direct interactions with the solvent, and (3) from the shift in spectral weight in the 0 K vibronic shape function. Although this second source of spectral broadening from NQEs is interesting, we cannot definitively say that the difference in the vibronic shape functions would hold if we were to compute them more rigorously.

### D. Comparison of E-ZTFC spectra with experimental spectra

Having computed both the vertical excitation energies of solvated snapshots and the ZTFC shape function, the absorption spectra including vibronic effects can be calculated in the E-ZTFC approach following Eq. (5). Combining our E-ZTFC approach with AI-PIMD sampling, Fig. 7, captures nearly all the experimental spectral widths along with the high energy tail.^{19,64} By contrast, using AIMD within the ensemble approach gives very poor agreement with the experimental spectrum with a full width at half maximum (FWHM) less than half that of the experiment (Table II). The spectrum computed with the E-ZTFC approach with AI-PIMD configurations doubles the FWHM. The inclusion of NQEs in the sampling causes 25% and 36% of this doubling of the FWHM for GFP and PYP respectively, with the remaining broadening arising from the inclusion of vibronic transitions using E-ZTFC.

. | GFP . | PYP . |
---|---|---|

Ensemble (AIMD) | 0.22 | 0.26 |

Ensemble (AI-PIMD) | 0.27 | 0.34 |

E-ZTFC (AIMD) | 0.37 | 0.40 |

E-ZTFC (AI-PIMD) | 0.42 | 0.48 |

Experiment | 0.51 | 0.64 |

. | GFP . | PYP . |
---|---|---|

Ensemble (AIMD) | 0.22 | 0.26 |

Ensemble (AI-PIMD) | 0.27 | 0.34 |

E-ZTFC (AIMD) | 0.37 | 0.40 |

E-ZTFC (AI-PIMD) | 0.42 | 0.48 |

Experiment | 0.51 | 0.64 |

For the remaining discrepancy between our computed spectra and the experimental spectra, one possibility is the error in the E-ZTFC approach. The use of an identical vibronic shape function for all vertical excitation energies in the ensemble is clearly an approximation within the E-ZTFC approach, as we found that Franck-Condon spectra varied between different frozen environments. It is difficult to quantify this error without computing a Franck-Condon spectrum for each snapshot in the ensemble and comparing the more rigorous resulting spectrum to the E-ZTFC approach, and even if we pursued this route, there would also be differences in the spectra due to the decreased sampling of anharmonic degrees of freedom of the chromophore in the rigorous approach.

Another possible reason for the remaining discrepancy between the E-ZTFC (AI-PIMD) results and the experimental data is that our snapshots are generated from a trajectory of the chromophore in pure water, whereas the experimental spectra are obtained in an aqueous solution of NaOH at pH = 13 for the GFP chromophore anion and a borate (III) buffer at pH = 10.2 for the PYP chromophore anion. None of the ions present in the experimental conditions were accounted for in our calculations because the size of the QM box used for the AIMD and AI-PIMD did not allow for realistic ion concentrations. The presence of the ions could provide additional electrostatic broadening to the ensemble spectra in these systems. An additional limitation of our MD simulations that may be responsible for some of the missing broadening is the small size of the AIMD and AI-PIMD periodic box, which could lead to solute-perturbed water structures near the periodic box edges.

Correctly capturing the breadth and tail of the experimental spectrum of the solvated PYP and GFP chromophores has been a challenge to previous theoretical approaches. For example, recent studies using implicit solvent and Gaussian broadening to account for inhomogeneous solvent effects^{72–74} required the use of a broadening parameter that was up to three times larger than that which would be predicted based on the solvent reorganization energy.^{73} This failure can likely partially be ascribed to the fact that the solvent reorganization energy was related to a solvent broadening parameter via the fully classical Marcus theory. It has previously been demonstrated that using a modified expression for the solvent broadening that accounts for NQEs yields a significant improvement in the computation of liquid phase photoemission spectra.^{75} Our results obtained for the PYP and GFP chromophore anions match this observation, as the AI-PIMD results in both cases predict a significantly larger effective solvent broadening than the AIMD data set. In our recent paper, in which we used a classical MM force field for the MD sampling,^{27} the FWHM of the simulated spectra for the GFP chromophore is 0.19 eV with the ensemble approach and 0.36 eV with the E-ZTFC approach, showing that AIMD does not lead to much additional spectral broadening compared to MM MD. By comparing the spectra computed with different approaches (see Sec. X in the supplementary material for further discussion), we have shown that the asymmetry and high-energy tail of the experimental spectra emerge from the high-energy excitation energies of direct chromophore-solvent intereactions, obtained from sampling the nuclear positions including NQEs, which are then further enhanced by the vibronic lineshape function.

## IV. CONCLUSIONS

In summary, we have performed an in-depth study of the influence of nuclear quantum effects, solvent effects, and vibronic effects on the computation of optical absorption spectra of two aqueously solvated chromophore anions. To analyze both direct and indirect effects of the solvent environment, AIMD and AI-PIMD simulations were performed for the chromophores in vacuum and with approximately one solvation shell of water around the chromophores. To include both short- and long-range solvent effects on the vertical excitations, TDDFT calculations on the MD configurations included one full solvent shell treated quantum mechanically and a large surrounding MM solvent region. The TDDFT vertical excitation energies computed from the ensemble of solute-solvent configurations can be efficiently combined with vibronic effects in our recently developed E-ZTFC approach that accounts for vibronic transitions via a zero-temperature Franck-Condon shape function computed in a frozen solvent environment. Absorption spectra were computed both in the ensemble approach and in the E-ZTFC approach.

We found that more accurate optical absorption spectra of the PYP and the GFP chromophores in water were computed if we explicitly account for the quantum nature of the nuclei, both in the AI-PIMD sampling of the ground state PES and through the inclusion of vibronic effects. The inclusion of NQEs in the sampling of the ground state PES through AI-PIMD yields significantly broadened vertical absorption spectra, as well as a spectral red-shift, compared to AIMD. From our analysis, we show that the influence of the AI-PIMD sampling on the computed optical spectra arises from three sources. First, in the chromophore alone, NQEs cause softening of heavy atom bonds whose motion couples strongly to the bright excited state. This softening leads to ensemble sampling from more anharmonic regions of the PES, causing a red-shift and spectral broadening. Second, AI-PIMD sampling yields an increase of solvent-induced spectral broadening for the vertical excitation energies, both indirectly through the configurations sampled by the chromophore and directly through solute-solvent interactions. However, although NQEs cause a strengthening of hydrogen bonds, our analysis indicates that the solvent-induced spectral broadening of vertical excitation energies is a collective effect that cannot be ascribed to specific hydrogen bonding sites. Finally, our results suggest that the AI-PIMD sampling may also influence the intensity of the vibronic transitions by leading to a redistribution of spectral weight from the 0-0 vertical transition to higher energy vibronic transitions. This redistribution of spectral weight leads to additional broadening of the absorption spectra.

Our analysis of the excitation energies performed on vacuum and solvated configurations shows that the presence of the solvent influenced the spectrum both indirectly, by affecting the configurations sampled by the chromophore, and directly through specific solute-solvent interactions. The spectral changes due to the indirect influence of solvent molecules suggest that both the chromophore and the explicit solvent should be fully flexible during the MD sampling of solute-solvent configurations.

Our computed spectral shapes show significant improvement in agreement with experimental spectra when using the E-ZTFC approach compared to the ensemble approach. The improvement is both in the spectral width and the high-energy tail. Although both including NQEs through sampling an ensemble of AI-PIMD configurations and including vibronic transitions through the ZTFC shape function leads to spectral broadening, the computed spectra remain slightly narrower than the experimental spectra for these aqueously solvated chromophores. We hypothesize that further increase in the computed spectral width may be achieved through a more rigorous coupling of solute and solvent in the calculation of the vibronic shape function, through the inclusion of counter ions in the MD or through the addition of homogeneous broadening due to the finite lifetime of the excited state of the chromophores in solution.

Overall, the E-ZTFC approach for computing optical spectra combined with AI-PIMD sampling of the ground state PES leads to a significant improvement in spectral shape compared to standard approaches. The combination of these methods is expected to yield much improved optical absorption spectra in a large variety of solvated dyes. The current work highlights the importance of NQEs when simulating optical absorption spectra, especially in systems where strong solute-solvent interactions are expected to occur.

## SUPPLEMENTARY METERIAL

See supplementary material for a detailed analysis of the nature of double-counting in the E-ZTFC approach, both in comparison to exact results for a series of model systems and in comparison to a double-counting free approach for the GFP chromophore. The supplementary material also contains additional convergence and validation tests regarding the computational protocol used to generate the solute-solvent configurations used in this work, as well as additional analysis regarding the structural changes in the chromophore due to the inclusion of NQEs and solvent effects. Furthermore, the structural changes responsible for observed variation in the vibronic shape functions for different frozen solvent conformations are described in detail.

## ACKNOWLEDGMENTS

This work was supported by the Department of Energy, Offce of Basic Energy Sciences CTC and CPIMS programs, under Award No. DE-SC0014437. T.E.M also acknowledges support from a Cottrell scholarship from the Research Corporation for Science Advancement and the Camille Dreyfus Teacher-Scholar Awards Program. This work used the XStream and MERCED computational resources supported by the NSF MRI Program (Grant Nos. ACI-1429830 and ACI-1429783). This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We would also like to thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that have contributed to these research results.

## REFERENCES

_{x}band of porphyrin as a case study