Modeling solvation effects on absorption and fluorescence spectra of indole in aqueous solution

The spectroscopic response of optically active molecules in a gas phase has been explored extensively over the years utilizing a broad spectrum of quantum mechanical methods that vary in accuracy and computational costs.^1–6 However, as most of the interesting photophysical and photochemical phenomena (i.e., photosynthesis, photocatalysis, UV absorption of DNA leading to skin cancer, and fluorescence in biological systems) occur in solution, in particular in aqueous environment, many efforts have been devoted to modeling the optical spectra of solutes in aqueous solution.^1,7–14 Ideally, one would like to compute excited states of the entire system (solute and solvent) quantum mechanically; however, due to the many-body problem that quantum mechanical methods entail, increasing the number of atoms in the system beyond a certain size results in an exceedingly large computational cost.¹³ Thankfully, each chromophore has a signature spectral response that is unique to the solute and remains present in the spectrum as the solute is embedded in different environments (i.e., gas and solution).¹ This feature allows one to focus on modeling the solvation effects on the optical spectrum of the solute rather than accurately modeling the solvent environment and its properties. Nonetheless, modeling the solvation effects on the optical spectra of molecules is still a challenge that computational chemists strive to overcome.

Features observed in the absorption and fluorescence spectra of a chromophore in solution include the position of the peak(s), especially that at λ_max, their corresponding intensities, and the broadness of the spectra. A large solvatochromic shift to the spectra is typically observed in solvents with high dielectrics such as water.¹ The peak intensities and shape of the spectra are also affected by the solvent. Some contributing factors are as follows: electrostatic and polarization effects, electronic density delocalization over the solvent molecules, and proton transfer reactions. Often, these factors are subtly combined when specific intermolecular interactions between the solute and the solvent, such as Hydrogen Bonding (HB), are present.¹³ Additionally, collective vibrational dynamics affect significantly the line shape.^12,15 From the modeling perspective, absorption and fluorescence spectra require that the solute ground state and excited state S₁ are calculated in the presence of the solvent.¹² For the sake of computational cost and efficiency, the solvent degrees of freedom are approximated in the theoretical models in a way that retains only the key components that contribute the most to the spectra. Thus, it is important to learn which of these effects are crucial to reproducing the various features of the absorption and fluorescence spectra and to identify a sufficient theoretical model that accurately captures these effects with minimal computational cost.

The Polarizable Continuum Model (PCM), an implicit solvation model, is a widely used approach to model the solvation effects on the electronic excited states inexpensively and at a high level of theory. In this approach, the solute is treated at a quantum mechanical level while immersed in a continuum polarizable dielectric medium mimicking the solvent.^12,13,16 While this approach can reproduce observed λ_max of the spectra in solution, it is not able to predict the spectral features, such as inhomogeneous spectral broadening, which result from explicit solute–solvent interactions.^1,12,17,18 A more explicit description of these interactions is obtained by using the hybrid Quantum Mechanics/Molecular Mechanics (QM/MM) approach or the QM cluster model. In the QM/MM approach, the system is partitioned into two regions: (i) the solute or QM region, treated at a high level of theory, and (ii) the solvent environment or MM region, treated at a lower but computationally efficient level of theory using force fields.^13,18 Using a non-polarizable force field, which is the most common approach, one can examine only some of the electrostatic effects due to the solvent on the electronic excited state of the solute. However, since force fields are based on empirical parameters fit to experimental observables or QM calculations, they may not be transferable enough to reproduce, for even moderately different systems, all the subtleties of solute–solvent interactions such as HB distances or the orientation of the first-shell waters.^1,19 To ensure an accurate prediction of the HB distance, orientation, and covalent contribution, a QM cluster model can be used where water molecules that are directly interacting with the solute are treated quantum mechanically. Using this approach also allows for mutual polarization between the solute and the solvent and has been found to improve on the prediction of the intensity at λ_max.²⁰ Nonetheless, it is computationally intensive, especially for systems with multiple HB sites. Moreover, it is limited to effects arising from short-range interactions since it ignores contribution from the bulk. Finally, like the PCM, both of these explicit solvation models are based on a static representation of the system and do not consider the dynamical nature of a solution.

At thermal equilibrium, both the solute and the solvent molecules are constantly in motion rather than occurring in a single static configuration. Sampling vibrational modes of the solute alone (in vacuum) using the Frank–Condon model couples the electronic and vibrational degrees of freedom, resulting in distinguishable transition bands in the electronic absorption spectra.¹² To capture the dynamical nature of the solvent molecules, it is customary to use molecular dynamics simulations; this additional layer of complexity is required to model the inhomogeneous spectral broadening observed in solution experiments. In particular, three classes of molecular dynamics simulations are typically considered: classical Molecular Dynamics (MD), in which nuclei are treated as classical particles interacting through classical potentials, ab initio Molecular Dynamics (AIMD), in which forces acting on classical nuclei are found by solving the Schrödinger equation, and Path Integral Molecular Dynamics (PIMD), in which nuclei are treated as quantum particles. Treating electrons or electrons and nuclei as quantum particles implies a progressively increasing computational cost, and thus, AIMD and PMID simulations are limited in terms of the size of the system and length trajectory; this limits severely the number of independent configurations sampled by the simulation.^12,13 In practice, an excited state calculation is performed for each uncorrelated configuration, and to simulate the absorption spectra in solution, the results are convoluted using a Gaussian or Lorentzian function. This method does not only account for the vibronic and dynamical effects but also includes thermal and pressure fluctuations that are also found to influence spectral broadening. Following this procedure, studies have found that modeling the vibrational effects and solvation dynamics is crucial to accurately reproducing this inhomogeneous broadening in the absorption spectra.^12,20,?,21 A more sophisticated approach includes vibrational effects of both the ground and the excited states using Franck–Condon overlaps.²² 

A crucial effect introduced by solvents, especially polar ones, is mutual polarization between the solute and the solvent molecules. Modeling this polarization effect in the ground state of a solvated molecule becomes computationally feasible with the development of polarizable force fields.^1,13,23 Amongst the many polarizable force fields that have been developed (i.e., induced dipole and Drude oscillator),²³ the most commonly used force field to study the mutual polarization effects on the absorption spectra in solution has been the fluctuating charge (FQ) model.^20,24–26 Modeling the polarization effects with this force field was found to successfully reproduce the absorption spectra of xanthine in solution with small deviation.²⁰ However, this was only achieved after fine tuning the atomic hardness parameter to exaggerate the polarization effects. Otherwise, the absorption spectra modeled with QM/FQ were found to be no different from that modeled with the non-polarizable force field, i.e., QM/TIP3P, as was also seen when studying the absorption spectra of curcumin tautomers.²⁶ Nonetheless, intensifying the polarization between the different solutes and water solvent indicated that it is important to model this effect in the ground state.

Capturing the mutual polarization effects becomes even more important in excited state calculations. For instance, upon UV-absorption, the chromophore becomes polarized by the light, and the electrons redistribute across the molecule. This change in the chromophore electronic density causes a rearrangement of the solvent charge distribution, which, in turn, polarizes the chromophore. Polarizable force fields combined with QM calculations in the QM/MM framework have been developed and used for this purpose.^14,27–31 Another approach that can include solvent polarization on excited states effectively and accurately is a hybrid approach using high level quantum mechanical methods combined with the Effective Fragment Potential (EFP) for the solvent.^19,32–35 In this model, the solute is treated quantum mechanically, while each solvent molecule is modeled as an effective fragment potential. In EFP, analytical potential functions such as the MM are used; however, the parameters are computed from first-principles methods and include multipoles beyond monopole (charges), including polarization (induction), dispersion, exchange and repulsion, and charge transfer terms in the interaction energy. Polarization with the excited states is included, which means that the solvent will be able to adapt to the change in the charge density of the solute upon excitation.¹⁹ 

One of the most studied chromophores, due to its biological relevance, is the optically active side chain of tryptophan, indole.^36–38 Indole is known for the sensitivity of its absorption and fluorescence spectra to changes in the microenvironment.^39–68 The lack of a detailed understanding of the molecular mechanisms resulting in its unique optical properties makes indole the perfect case study for investigating the solvation effects on the absorption and fluorescence spectra. We focus on investigating how the different solvation effects modeled by PCM, QM cluster, and explicit solvation models contribute to the accuracy of the predicted absorption and fluorescence spectra. In particular, since explicit solvent models entail sampling of molecular configurations of the electronic ground state, we adopt the molecular dynamics simulation protocol using force-field or Density Functional Theory (DFT)-based Car–Parrinello schemes to integrate the equations of motion. Both non-polarizable and polarizable (AMOEBA) force fields are considered, while for Car–Parrinello molecular dynamics, we use a plane wave expansion of the Khon–Sham orbitals and norm-conserving pseudopotentials. The excited states of indole in the various solvation models are studied with the equation of motion coupled cluster with single and double excitations (EOM-CCSD) and the EFP method.

Previous investigations using explicit solvent models have focused on absorption spectra. Studies on fluorescence exist as well;^69–73 however, modeling fluorescence spectra is challenging because calculating the gradients of the excited states’ potential energy surfaces can be complicated. Additional problems occur when more than one state absorbs light, as is the case in indole, and predicting the correct emitting state becomes an important issue. Here, we offer an initial investigation into these issues in indole. To explore the configurational space of an excited indole moiety in aqueous environment, we develop novel force-field parameters to be used with the SPC water model and calculate the solvent reorganization free energy following chromophore excitation. We also explore the capability of our models to predict the correct emitting state. The current study offers an insight into additional further steps that are needed for modeling fluorescence.

Geometry optimization and vibrational frequency calculations were performed on indole in vacuum at the B3LYP/6-31G(d) level of theory. Four vertical excitation energies were computed from the ground state minimum with EOM-CCSD/6-311+G(d) in aqueous solution implicitly using the PCM with the Perturbation Theory Energies and Singles (PTES) scheme.^74–76 Due to the fast time-scale that absorption occurs (≈10⁻¹⁵ s), the nuclei of the solvent do not have time to respond to the change in the electronic density of the solute. Instead, the fast electrons of the solvent, represented by the dielectric constant, respond to the polarized excited state density from which the excitation energies are computed using the non-equilibrium solvation approach. The excited state relaxation to the S₁ minimum of indole in aqueous solution was then computed by optimizing with EOM-CCSD/6-311+G(d) and PCM-PTES with a state-specific approach. Since fluorescence happens at a slower timescale than absorption (≈10⁻⁹ s), the slow nuclei of the solvent have time to reorient and align their dipoles to those of the solute in the excited state, a process also known as solvent reorganization. Here, the solvent becomes equilibrated with the electronic density of the specified S₁ state of the solute, and the energy is resolved using the equilibrium solvation approach.

To evaluate how HB influences the excitation energies and intensities of absorption, an indole +1H₂O complex was studied. We modeled this complex by placing 1H₂O at the N–H site of indole where it would be expected to hydrogen bond (as shown in Fig. 1).

The complex was optimized to the ground state minimum at the B3LYP/6-31G(d) level with the addition of Grimme dispersion correction, D3, to account for the dispersion error from TDDFT methods, which improves the prediction of the HB interaction. A vibrational frequency calculation was performed to confirm that it was a true local minimum. Four vertical excitation energies from the equilibrium complex in the ground state were computed at the EOM-CCSD/6-311+G(d) level to study the absorption. Excited state optimization and vibrational frequency calculations with CAM-B3LYP+D3/6-311+G(d) were performed to obtain the S₁ minimum from which fluorescence occurs. All geometry optimization and vibrational frequency calculations were performed on the Gaussian16 software.⁷⁷ A single point calculation was performed at the optimized S₁ minimum of the indole +1H₂O complex using EOM-CCSD/6-311+G(d) to obtain the vertical emission energy and corresponding oscillator strength.

The vertical excitation and emission calculations were then repeated using the hybrid Quantum Mechanics and Effective Fragment Potential (QM/EFP) approach,^19,33–35 with the optimized ground and excited state geometries of the indole +1H₂O complex. Indole was described as the QM region computed at the EOM-CCSD level, and the 1H₂O molecule was modeled with EFP with parameters computed at the same level of theory.

A structural model of indole solvated by a cubic box of 3291 water molecules was simulated using Gromacs 2016.3.⁷⁸ The equilibrium geometry for indole was calculated from energy minimization in the gas phase at the B3LYP/6-31G(d) level of theory. Fixed atomic charges were computed using the ChelpG approach⁷⁹ on the ground state equilibrium geometry at the EOM-CCSD level of theory. All the other force field parameters, including that of the 12-6 Lennard Jones (LJ) potential used to model the van der Waals (vdW) interactions, were taken from the Generalized Amber Force Field (GAFF) using the AmberTools Leap program.⁸⁰ We used the SPC model for water molecules.⁸¹ The benchmark of the atomic charges can be found in the electronic supplementary material (ESI).

The geometry of the system was initially optimized using the steepest-decent algorithm before starting the molecular dynamics simulation. The simulation included Periodic Boundary Conditions (PBCs) along the xyz directions. The smooth Particle Mesh Ewald (PME) method was used to compute long range electrostatic force terms with a grid spacing of 0.16 nm. This consisted of an equilibration phase (lasting for 10 ns) at constant volume in which the system was allowed to reach 300 K using the V-rescale, modified Berendsen thermostat (the NVT ensemble) with a temperature coupling tau t of 0.1 ps.⁸² The volume was then allowed to fluctuate by setting the pressure to 1 atm in two steps using first the Berendsen and then the Rahman–Parrinello barostat algorithms^83,84 with a pressure coupling tau p of 2.0 ps. After 10 ns, the system was equilibrated, and thus, we started the MD production phase (hereafter referred to C-MM/MD simulation) during which we sampled the ground state configuration space for a total of 10 ns. In equilibration and during the C-MM/MD steps, the length of all the bonds was kept fixed using the LINCS algorithm,⁸⁵ and the equation of motion was integrated using the leapfrog⁸⁶ using a time step of 2 fs. A subsequent run of 5 ns was performed by releasing the constraint on the bond length and decreasing the time step to 0.5 fs. This is denoted as UC-MM/MD.

To evaluate the sampling timescale needed to reproduce structural features of the solvated system, 20 uncorrelated snapshots separated by 250 ps were extracted from the 5 ns UC-MM/MD simulation and used as the random initial guesses to run short 8.3 ps MD simulations with a time step of 0.05 fs. The first 5.3 ps of the simulations were considered equilibration and the remaining 3 ps as the MD production used for analysis.

A cubic box including one indole molecule and 3991 waters was modeled using the polarizable AMOEBA force field and the Tinker 8.7 molecular dynamics package.⁸⁷ Parameters from amoeba09⁸⁸ and amoeba03⁸⁹ were utilized for the potential energy terms that make up the force field for indole and the water molecules, respectively. Unlike the non-polarizable force field, AMOEBA uses the 14-7 pairwise additive function for the vdW interaction and includes both permanent and induced multipole moments at each atomic center to calculate the electrostatic interaction.⁹⁰ In addition, electronic polarization is introduced in the force field through the mutual induction energy term, which allows each polarizable atomic site to respond to an external electric field. The induced dipole at any atomic center will further polarize other sites until all atoms in the molecule have been mutually polarized.^90,91 The large network of water molecules act as an electric field around the indole molecule. Modeling both indole and the water molecules with AMOEBA, thus, allows indole to become polarized by the surrounding waters and vice versa.

The system was optimized using the steepest-decent minimization algorithm and then equilibrated. The PME method was used to compute electrostatic interactions with a real space cutoff of 0.70 nm and periodic boundary conditions along xyz. The equilibration occurred in two steps using the NVT and NPT ensembles. The system was heated to 300 K with the Berendsen thermostat and a tau of 0.1 for 50 ps. The pressure of the system was then equilibrated to 1 atm by simulating the system for 5 ns with a dt of 1.0 fs. This was done using the Berendsen barostat using a pressure coupling tau p of 1.0 ps. The equilibrated system was then used in the MD production phase and further simulated for 5 ns with a dt of 0.5 fs using the velocity Verlet integration method.

A cubic box of one indole and 126 waters was modeled with a plane wave and pseudopotential and Strongly Constrained and Appropriately Normed (SCAN) functional of Density Functional Theory (DFT) on quantum espresso 6.4.1.⁹² Norm-conserving HSCV-PBE pseudopotentials⁹³ were used for the core electrons and nuclei of all atoms involved (N, C, O, H). A kinetic energy plane wave cutoff of 85 Rydberg was used (the benchmark can be found in Fig. S5 of the supplementary material). A fictitious electronic mass of 300 a.u. was used. The system was then simulated for 8.3 ps with a dt of 0.05 fs using Car–Parrinello Molecular Dynamics (CPMD), of which 5.3 ps were considered equilibration and the last 3 ps were considered production and used for analysis.

Four excited states were computed for 50 snapshots of indole with water molecules within 4 Å cutoff distance sampled from the C-MM/MD, UC-MM/MD, AMOEBA MM/MD, and SCAN-CPMD simulations, illustrated in Fig. 2(a). The solvation shell with a cutoff of 4 Å included an average of 18 ± 2 waters for configurations modeled with force field simulations and 15 ± 1 water molecules for those modeled using SCAN-CPMD. The 50 snapshots were separated by 200 ps, 50 ps, and 50 ps and 40 fs for the sampling methods, respectively. To evaluate the effect that long-range solute–solvent interactions have on the absorption and check the accuracy in describing the unique HB network between water molecules in the liquid phase, four excited states were computed for the same 50 snapshots from the UC-MM/MD simulation including 6 waters and 8 Å of indole, as shown in Figs. 2(b) and 2(c). Discussions on the effects of different cutoff radii are given in the supplementary material and show that although there are some differences when increasing the radius and the number of water molecules, the main picture is captured with the small radius.

All excited state calculations of indole and water clusters sampled from MD simulations were modeled using the EOM-CCSD/EFP approach with basis set 6-311+G(d) for high accuracy at an affordable computational cost. The EFP model allows for an efficient and accurate description of the waters surrounding the solute and is described in detail elsewhere.^19,32–35 In the coupled-cluster calculation with EFP, the induced dipoles of the fragments are frozen at their HF values, so mutual polarization is introduced with the ground state. To account for solvent’s response to electron rearrangement in the excited states, a perturbative non-iterative correction is computed for each EOM root. The one-electron density of the target EOM state is calculated and used to re-polarize the environment, i.e., to recalculate the induced dipoles of the EFP part in the field of an EOM state. These dipoles are used to compute the polarization energy corresponding to this state, and this polarization energy is added to the energy of the excited state.

The calculations were repeated with CAM-B3LYP+D3/6-311+Gd(d) for models with a centered indole solute surrounded by waters within a 4 Å radius in order to test the effect of a different electronic structure method. All excited state related calculations were computed on QChem.⁹⁴ The absorption spectra were simulated from the four excitation energies and corresponding oscillator strengths with a Lorentzian function with a bandwidth of 0.2 eV.

In order to validate the EFP approach, we calculated the spectra using complexes of indole with a single water molecule either all quantum mechanically or using EFP for the water. The comparisons are given in the supplementary material and show comparable results, validating our approach of using EFP for the larger clusters.

We modeled indole in the excited L_b and L_a states using the simple non-polarizable force fields with classical molecular dynamics in order to study the fluorescence of indole in an explicit solvation model. A proper treatment of using MD for the excited states requires reparametrization, and there have been several efforts to develop force fields for excited states.^95–100 Here, we used the C-MM/MD model that is used for sampling the ground state and parameterized the non-bonding force field potential terms, the vdW potential and the Coulombic potential for the excited L_b and L_a states of indole, respectively.

The Coulombic potential was parameterized by computing the partial charges of an isolated indole molecule using the densities of L_a and L_b in the ground state equilibrium geometry with ChelpG at the EOM-CCSD/6-311+G(d) level of theory. Ideally, to properly model the fluorescence, these charges should be computed at the respective excited state potential energy surface minima; however for indole, the energy gap between the L_b and L_a states at emission is too small, leading to mixing between the two wavefunctions. This generates charges of the two excited states, which do not reflect their dramatic difference in polarity. This mixing does not occur for the geometry corresponding to the ground state potential energy surface minimum, and that is the reason why we use this geometry to calculate charges for the excited states L_b and L_a. The derived charges can be found in the supplementary material (Table S1).

Non-polarizable force fields use the 12-6 LJ potential to describe the van der Waals interactions,

where r is the distance between particles of type i and j, respectively. For each pair of particle types (i, j), this potential is completely determined by parameters ϵ_i,j and σ_i,j or, equivalently, by pairs C₁₂(i, j) and C₆(i, j). Here, we consider how the electronic excitation affects these parameters and, in particular, C₆(i, j), which according to the Slater–Kirkwood equation¹⁰¹ depends on atomic polarizability as follows:

where E₀ is the ionization energy, a₀ is the atomic radius, α_i is the polarizability (in units of V⁻¹), N_A is the Avogadro number, and N_eff,i is the number of effective electrons in atom i. Since the atomic polarizability α_i is directly proportional to the atomic volume,²¹ the denominator in the right-hand side of Eq. (2) is proportional to $ViNeff,i$⁠, where V_i is the atomic volume, and it is, therefore, also proportional to $1ρi$⁠, where ρ_i is the atomic density,

The effect of the electronic excitation is, thus, to rescale the attractive C₆(i) pair of the ground state,

by a factor k equal to

where $C6s$⁠, V_s, and ρ_s are the C₆ coefficient, the atomic volume, and the atomic density for the electronic state s, respectively. We computed the charge density of the atoms in indole for the ground state and the two excited states at the EOM-CCSD/6-311+G(d) level. Atomic charges and atomic volumes were then computed using the Bader analysis grid-based approach.¹⁰² The scaling factor k was calculated for each atom in indole. The resulting ϵ parameters used in the MD simulation for each atom in the ground and excited states are reported in the supplementary material (Table S2).

In the gas phase, indole is found to be emitted from the L_b excited state. Upon solvation, the L_a excited state is redshifted below the L_b excited state, becoming the S₁ minimum where fluorescence takes place. We model this process (as shown in Fig. 3) by performing the optimization and vibrational frequency calculation of the electronic ground and L_b and L_a excited state minima of indole in the gas phase, quantum mechanically, at the TDDFT level. A thermal correction to the Gibbs free energy is obtained from the vibrational frequencies and added to the EOM-CCSD electronic energy to obtain the free energy of the solute in the respective electronic states. The difference between the free energy of the solute in the ground and excited state minima in the gas phase is defined as the adiabatic free energy in the gas phase, $ΔGgasadiab$ (depicted by a blue solid arrow in Fig. 3).

To evaluate whether modeling the L_b and L_a excited states of indole in explicit solvation with the L_b-MM/MD model and L_a-MM/MD models can predict this energy level inversion, we calculate the adiabatic free energy of the two excited states in solution, $ΔGsoladiab$ (depicted by the blue dotted arrow). We estimate this free energy difference as resulting from three steps: (i) desolvation of indole in the ground state, (ii) excitation of indole in the gas phase and relaxation of the indole’s nuclear degrees of freedom (this allows us to estimate $ΔGgasadiab$⁠), and (iii) solvation of indole in the excited state (Fig. 3).

The adiabatic free energy in solution, $ΔGsoladiab$⁠, is also informative about solvent reorganization free energy. Upon excited state relaxation to the excited state minimum in solution, solvent molecules reorganize in response to the change in the charge distribution of the solute. The free energy required for solvent molecules to reorient in response to a change in the solute charge distribution (as a result of electronic excitation) is known as solvent reorganization free energy, ΔG^reorg. We define ΔG^reorg as follows:

We compute the free energy of solvation, ΔG^sol$v$, starting with the (GS) C-MM/MD, L_b-MM/MD, and L_a-MM/MD models of indole in the ground and excited states in aqueous solution using the Free Energy Perturbation (FEP) theory and molecular dynamics approach. We decouple indole from water by decreasing the value of a coupling parameter, λ, that modulates the strength of the Hamiltonian terms $ΔĤ$ associated with the interactions between the solute and the solvent. We used a decoupling protocol in which λ is varied from 1 to 0 in decrements of 0.05 for vdW terms and by decrements of 0.1 for the Coulombic terms. Decoupling of electrostatic terms starts after vdW interactions have been scaled by a factor of 0.55. The FEP/MD simulation was run for 100 ns for each λ value to enable statistical convergence of the free energy estimated by combining the 20 simulations (each with a different value of λ) (see the supplementary material).

From the L_a-MM/MD simulation, 50 snapshots separated by 200 ps were extracted. The frames included water molecules within 4 Å radius of indole, and the indole geometry in all 50 configurations was superimposed and replaced with the quantum mechanically optimized L_a geometry. Single point excitation energies were computed for all 50 frames with EOM-CCSD/EFP and 6-311+G(d) basis sets. The fluorescence spectra were simulated from the computed S₁ emission energy and corresponding oscillator strength dressed with a Lorentzian function with a bandwidth of 0.2 eV.

The character of all excited states computed during absorption and emission with all models described was characterized from Natural Transition Orbital (NTO) calculations.

We begin by investigating the structural features of the indole electronic ground state in aqueous solution. We focus on the radial distribution function g_H–O(r) calculated with UC-MM/MD, AMOEBA MM/MD, and SCAN-CPMD simulations. These methods show a hierarchy of approximations in terms of electrostatic interactions: from fixed-charge particles (UC-MM) to point dipoles (AMOEBA) to a full description of the electronic charge density (SCAN). As an additional term of comparison, we also consider the results from an intermediate level of approximation based on the BLYP exchange and correlation functional and Born–Oppenheimer MD. BLYP is based on a Generalized Gradient Approximation (GGA), which is, in general, less accurate than the meta-GGA scheme utilized in SCAN. Data for the radial distribution function using BLYP Born–Oppenheimer MD are taken form Ref. 103. Since an experimental g_H–O(r) is not available, we cannot characterize the accuracy of each method. Instead, we discuss similarities and differences between different theoretical models. g_H–O(r) are plotted in Fig. 4.

The peak in g_H–O(r) obtained with UC-MM/MD indicates a typical distance of 1.875 Å between the amino hydrogen of indole and the water oxygen. This is consistent with the optimized geometry of the indole +1H₂O complex for which the distance is found to be 1.88 Å at the CAM-B3LYP+D3/6-311+G(d) level of theory. On the other hand, both the AMOEBA MM/MD and SCAN-CPMD simulations predict the water oxygen at a larger distance from the amino hydrogen of indole, specifically at r = 2.025 Å, in agreement with the results of the BLYP-BOMD simulation. This discrepancy between the optimal distance found for the indole +1H₂O complex and the peak of the g_H–O(r) calculated with SCAN likely results from competitive interactions due to the network of H-bonded water molecules from the second shell, which “pull” the first shell water away from the minimum of the potential energy surface. This exquisite many-body effect may not be present in the classical non-polarizable force field whose parameters are typically optimized to reproduce the geometry of a supramolecular complex with a water molecule. Interestingly, the inclusion of atomic polarizabilities in AMOEBA⁹⁰ results in a more accurate structure with a typical distance between amino hydrogen and water oxygen completely consistent with the DFT-MD calculations.

The three simulations (UC-MM/MD, AMOEBA-MM/MD, and SCAN-CPMD) also disagree on the predicted height of the first peak and position and magnitude of the first minimum. For UC-MM/MD and AMOEBA MM/MD, g_H–O(r) shows peak values of 1.226 and 1.026 and decays to 0.250 and 0.413 at a distance of 2.475 Å and 2.625 Å, respectively. g_H–O(r) predicted with SCAN-CPMD shows an even more pronounced peak (1.368) and decays to a minimum of 0.07 at 2.475 Å. Overall, SCAN-CPMD shows a significantly structured first solvation shell compared to the other methods. Despite the stark contrast between SCAN-CPMD and the other two simulations, we refrain from drawing conclusions, given the significantly different timescales (3 ps in SCAN-CPMD vs 5 ns for the force-field based ones) sampled along the trajectory. Integration of the radial distribution function over the first peak indicates a coordination number of ∼1 for all the simulations. This means that all methods agree on the fact that only one water molecule is HB to indole at the NH site. This further justifies our direct comparison of the HB distance in the indole +1H₂O complex with that predicted from the MD simulations.

The UC and AMOEBA MM/MD simulations do not show pronounced second peaks in g_H–O(r), indicating a substantial lack of the structure in the second shell of solvation around the NH group of indole. By contrast, the SCAN-CPMD simulation predicts a peak at about 3.75 Å with a g_H–O(r) of 1.22 (almost as large as the first peak) with the second local minimum located at about 5 Å, indicative of the presence of a significantly structured second solvation shell. g_H–O(r) predicted by the BLYP-BOMD simulation shows a somewhat intermediate character with a broad second peak located between 2.8 Å and 6 Å. Despite the fact that the SCAN-CPMD simulation is based on a relatively accurate description of the electronic degrees of freedom and that computed g_O–O(r) (see Fig. S6 of the supplementary material ) gives best agreement with the experimentally measured liquid water structure, the length of the simulation is too short to reach the statistical convergence of estimated g_H–O(r), and as a result, some of features of g_H–O(r) apparent from the plots may result from statistical fluctuations. However, some structural differences appear to be robust and will be further explored below in terms of their consequences on the predicted absorption spectra.

In order to check the role of the simulation time, we ran 20 classical simulations that are the same timescale and time step as that of the Car–Parrinello simulation starting from 20 uncorrelated snapshots extracted from the 5 ns UC-MM/MD trajectory. We plotted the radial distribution function of the statistical average and standard deviation at each separation distance r and compared the results to the RDF of the 5 ns UC-MM/MD simulation and the 3 ps SCAN-CPMD simulation (see Fig. S7 of the supplementary material). Our results indicate that the difference in the first and second peaks between classical and CPMD simulations is unlikely to be due to statistical fluctuation. It also indicates that the differences we see in the spectra in SCAN-CPMD compared to the UC-MM/MD simulation are due to the short length of the simulation. These results also indicate that a CPMD study could be devised that could give statistically converged results. This is, however, computationally very intensive, prohibiting the application to the current study.

Absorption spectra were simulated using four vertical excitation energies and oscillator strengths of 50 indole–water frames sampled with the non-polarizable force field with constrained (C) and unconstrained (UC) bonds, the polarizable AMOEBA force field, and SCAN-CPMD simulation. Waters in the solvation shell around indole were selected using a 4 Å cutoff radius, and vertical excitation energies were computed with the EOM-CCSD/EFP approach. For comparison, the absorption spectrum was also simulated using the optimized configuration of indole with the EOM-CCSD/PCM-PTES approach. All spectra were overlaid with the experimental spectra of indole in aqueous solution plotted in Fig. 5.

The experimental spectrum exhibits two peaks, at 270 nm (4.59 eV) and 285 nm (4.35 eV), with the peak at higher energy being broad and having some underlying structure. We will investigate how the different solvation models perform in accurately reproducing the main features of the absorption spectrum of indole in aqueous solution: the position of λ_max at 4.59 eV, the position of the second peak at 4.35 eV and its relative intensity, and the spectrum line shape.

Since the focus of this work is on the solvation effects in order to compare the theoretical spectra to the experimental one, we have shifted them by 0.46 eV. This corresponds to the error of the EOM-CCSD/6-311+G(d) method to reproduce the gas phase vertical excitation energy of indole according to our previous work.⁴³ We first examine how the different solvation models predict the position of the λ_max peak. Modeling the solvation effect on the absorption of indole implicitly with the PCM is found to predict λ_max at 4.83 eV, blueshifted by 0.24 eV compared to the experiment. It should be noted that since we have already extracted the EOM-CCSD error, this shift and all other ones discussed in this section are purely an error of the solvation description. Explicit solvation while keeping the solute structure rigid in the C-MM/MD model gives almost the same λ_max, while the peak becomes broader. There is no surprise that the PCM-PTES and C-MM/MD models predict similar λ_max since they both use the same frozen optimized ground state geometry of indole. The broadness of the C-MM/MD peak reflects the effect of explicit solvent interactions with indole. Including vibrational effects in both indole and the water molecules in UC-MM/MD by relaxing the bonds also does not change λ_max by more than 0.03 eV. This demonstrates that the vibrational effects via classical molecular dynamics (excluding any QM vibronic effects) have little influence on the position of λ_max. The other two models, however, which include mutual polarization between all the molecules, produce greater shifts, an indication that the polarization effects in the ground state are crucial to reproducing the position of λ_max. SCAN-CPMD predicts a blueshift in the position of λ_max by 0.08 eV compared to the experiment, whereas AMOEBA MM/MD finds it to be redshifted by 0.17 eV from the experimental peak.

Both the experiment and most theoretical spectra predict that the first peak has a shoulder, although the relative intensities and positions differ. We will discuss what leads to the λ_max and local maximum peaks and what leads to the differences predicted by the MD simulations in more detail in Sec. III B 4.

We can also compare to the gas phase λ_max and extract the solvatochromic shift for this peak from the different models and compare to the experimental value, which is redshifted by 0.28 eV.^42,43 All methods predict a redshift, and the magnitude ranges from 0.1 eV to 0.4 eV. The PCM gives the smallest value of 0.1 eV, while AMOEBA predicts the largest shift, about 0.4 eV. The other methods predict a shift of about 0.2 eV. Here, we have estimated the solvatochromic shift in our calculations by comparing λ_max to the vertical excitation energy in the gas phase. This estimate assumes that λ_max in the gas phase corresponds to the vertical excitation energy.

In the measured absorption spectra, a second peak is observed at 4.35 eV with a relative intensity of 0.72 compared to the first peak, and as it is not well defined, it still appears as a pronounced shoulder in the spectrum. Both PCM and C-MM/MD predict the second peak to be at about the same position (4.42 eV and 4.44 eV, respectively) for the same reason as in their similarity in λ_max. They, however, differ in the relative intensity. The relative intensity of the second peak is predicted to be 0.29 by the PCM. C-MM/MD enhances this intensity to 0.86. Including the vibrational effects empirically from the UC-MM/MD model somewhat redshifts the peak, while it increases its intensity. SCAN-CPMD portrays this peak to be at 4.36 eV, which is only 0.01 eV away from the experimental value; however, the relative intensity is underestimated. AMOEBA predicts the second peak to be at 4.28 eV, redshifted by 0.07 eV compared to the experiment. In general, force field approaches overestimate the intensity of the second peak, while SCAN-CPMD and PCM underestimate it. Overall, the position of the two peaks is best described by SCAN-CPMD, but the structure and overall shape are better described by AMOEBA.

In order to better examine the shape of the spectra, in Fig. 5(b), computed λ_max of the simulated absorption spectra were all shifted to the experimentally observed peak at 4.59 eV. Although SCAN-CPMD accurately reproduces the position of the two peaks [λ_max and the second peak, as shown in Fig. 5(a)], it fails to reproduce the observed spectrum broadness, sharp decay from the second peak, and the low intensity–low energy tail trail. SCAN-CPMD portrays a spectrum that resembles the shape of that produced from the PCM and C-MM/MD simulation, both of which do not include any vibrational flexibility. The spectrum is broader than that in the PCM yet narrow in comparison to the C-MM/MD and the experimental spectrum. It should be noted that a snapshot from the C-MM/MD simulation was used to start the simulations using the SCAN-CPMD approach. The similarity in the spectral shape of SCAN-CPMD to the PCM and C-MM/MD suggests that perhaps running the CPMD simulation for 5.3 ps to equilibrate is not enough to decorrelate the system from the reference state and/or that snapshots that are 40 fs apart are too correlated, more uncorrelated configurations need to be sampled. This confirms the observation and discussion we had in Sec. III A 2. On the other hand, both non-polarizable and polarizable force fields seem to reproduce the broadness of the spectra, although they overestimate it because of overestimating the low energy peak. AMOEBA is the model that best reproduces the shape of the experimental spectrum, including the width and the shape at both low and high energy ends of the spectrum.

We explore the underlying electronic transitions responsible for the predicted features seen in the absorption spectrum of indole in solution. The absorption spectrum of indole in this energy range is governed by excitation to two ππ^* excited states, called L_a and L_b. Two πσ^* states also exist nearby. In order to find the contribution of each state to the spectrum, we strip down the absorption spectra from the UC-MM/MD, AMOEBA MM/MD, and SCAN-CPMD simulations into three spectra simulated from the energies and oscillator strengths corresponding to the electronic L_b, L_a, and two πσ^* state transitions. Figures 6(a)–6(c) show the individual spectra superimposed to the total spectrum in each case. Furthermore, in order to compare the effects due to vibrations from the effects of interactions with water molecules, we repeat the calculations and plot the spectra using indole only without the waters [Figs. 6(d)–6(f)]. It should be noted, however, that these are not gas phase spectra since they are obtained from the condense phase simulations, so the vibrations have been affected by the surrounding solvent molecules.

It is obvious from Fig. 6 that the two main peaks seen in all spectra arise from the electronic states L_a and L_b. The higher energy peak that corresponds to λ_max is due to the L_a state, while the second peak at lower energy is due to the L_b state (mostly). It is also clear that while the L_b band has a single maximum and is described well by a Lorentzian, this is not the case for the L_a band. It shows a double maximum or a shoulder, and as was discussed above, these features are apparent in the overall spectrum. In the UC-MM spectrum, the lower energy band is actually a combination of both the L_a and the L_b states since one of the subbands in L_a blends with the L_b band. The combined band is broader, and the maximum is shifted to higher energies compared to the L_b maximum. This highlights the complexity of assigning the bands of the absorption spectra and the potential errors when we compare only to vertical excitations from calculations.

The πσ^* states are also present at energies near L_a, and they likely contribute to the overall intensity and shape of the peaks. In order to further probe their contributions, we plotted the spectra with and without them (shown in Fig. 7). Figure 7 confirms that while the overall spectra do not change, the relative intensities are affected by the πσ^* states, primarily at the high energy end of the spectrum. This is particularly true for UC-MM, where the peak with the higher intensity is no longer the one at the high energy when the πσ^* states are excluded. The effect on the other two simulation methods is much smaller.

In order to explore the origin of the detailed structure in the L_a band, we examine the vibrational and solvent interaction effects separately. We examine how the vibrational modes of indole contribute to the absorption spectra by computing the excited states using only indole without the surrounding water molecules [Figs. 6(d)–6(f)]. In the gas phase, the L_a spectra show a shoulder in the main peak in the UC-MM/MD and AMOEBA MM/MD simulation. In the UC-MM spectra, the shoulder is separated by 1000 cm⁻¹ from the peak [panel (a) of Fig. 6]. In the L_a spectra simulated from the AMOEBA MM/MD configurations [panel (b) of Fig. 6], the separation is 1033 cm⁻¹. SCAN-CPMD, on the other hand, only predicts a single peak [panel (c) of Fig. 6]. It is clear from these spectra that the L_a band exhibits a vibrational progression, while L_b does not. This is in agreement with their minimum energy geometries and how they differ from the ground state minimum. The L_b minimum is very similar to the ground state minimum, so the Franck–Condon factors do not favor a vibrational progression. The L_a minimum, on the other hand, differs significantly from the ground state minimum along the C–C and C–N bonds, and this will lead to vibrational progressions based on Franck–Condon arguments. The reason that SCAN-CPMD does not display any vibration is likely because there was not enough time for vibrational relaxation.

We now explain how the interaction with the water environment influences the spectra produced from each simulation [Figs. 6(a)–6(c)]. Addition of the water interactions mainly changes the relative intensities of the L_a and L_b peaks. It also increases the separation of the vibrational progressions in L_a. In UC-MM, the separation between the two peaks is 1722 cm⁻¹, while in AMOEBA, it is 1301 cm⁻¹. In AMOEBA, the gap between the L_a and L_b peaks decreases with solvation to the extent that L_b becomes almost completely hidden beneath the L_a spectra [panel (e) of Fig. 6].

The πσ^* states have a much more pronounced effect with the addition of water molecules. Two distinct peaks for the πσ^* spectra become less structured and form many peaks, an indication that the energy of the πσ^* states is very sensitive to the position and orientation of the surrounding water molecules, as one would expect.

In order to examine the effect of different electronic structure methods, the spectra were also calculated using CAM-B3LYP+D3/6-311+G(d) using the same configurations. The simulated absorption spectra are shown in Fig. 8. These spectra are not shifted, so it is obvious that they reproduce the positions better than EOM-CCSD. On the other hand, the structural features are diminished. The two distinct peaks simulated from the various MD models and the PCM at the EOM-CCSD level seen in Fig. 5 are less obvious in these spectra, and sometimes a second peak completely disappears. This is because the energy gap between the two states is, in general, underestimated by the CAM-B3LYP method, so there is no distinct peak separation. The position of λ_max and spectral broadness of the absorption spectra simulated from the AMOEBA MM/MD approach best agree with the experiment, relative to other ground state sampling methods. This again points to the importance of modeling the mutual polarization effects in the ground state and the capability of polarizable empirical force fields such as AMOEBA to properly sample solute–solvent interactions. Although the spectra predicted at the TDDFT level may appear to better predict the position of the spectra, the underlying transitions that make up the features of the spectra are incorrect. Characterizing the excited states from NTO calculations portrays that simulating the absorption of indole with an implicit solvation model, such as the PCM at the CAM-B3LYP level of theory, predicts the L_b(1ππ^*) state to be brighter than L_a(2ππ^*). Modeling the solvent explicitly, on the other hand, predicts that indole absorbs mostly to the S₁ excited state and not S₂. It also predicts the wrong ordering of states at the ground state minimum, with the S₁ state being L_a and the S₂ excited state being L_b. This demonstrates that while CAM-B3LYP is a favorite TDDFT functional for computing excited states of neutral systems, it is inadequate for modeling the solvation effects both implicitly and explicitly for indole.

The performance of other TDDFT functionals and basis sets in accurately predicting the excited states of indole in explicit solvation is not tested since that is not the focus of the current paper. Several TDDFT functionals have been tested for indole in the gas phase and aqueous solution using the implicit solvation model called the PCM.³⁹ A benchmark study with various TDDFT functionals and basis sets would be necessary to identify or confirm whether or not TDDFT is a reliable method for indole in solution (explicitly) more generally.

Upon excited state relaxation to the S₁ minimum where fluorescence takes place, the water molecules rearrange from their initial orientations at the ground state minimum so that their dipoles align optimally in response to the change in electronic and nuclear charge distribution of the solute in the excited state minimum. For indole, the polar L_a state is found to be largely stabilized by the water solvent to the extent that an energy inversion between the L_a and L_b excited state takes place at the fluorescence where L_a becomes the emitting state.^{40,48,104–106} Hence, experimentally, it is observed that while in the gas phase the emission occurs from the L_b state, the ordering of the states is inverted in aqueous solution and emission occurs from L_a. We will examine the ability of implicit and explicit solvation models to predict this level inversion and accurately reproduce the fluorescence spectra of indole in aqueous solution.

We model indole in the excited state in aqueous solution using the explicit non-polarizable SPC water solvent and molecular dynamics (similar to the C-MM/MD model used to model indole in the ground state for absorption). We modeled the excited L_b and L_a states of indole by parameterizing the non-bonding parameters of the force field, as described in Sec. II. Water molecules are equilibrated to the different electrostatics of indole in the respective electronic states using MD simulations, which we term L_b-MM/MD and L_a-MM/MD. The changes in the solvation structure upon excitation are shown by plotting g_H–O(r) for the ground and excited L_b and L_a states (shown in Fig. 9).

The HB strength with water becomes stronger in the L_a state, as is indicated in the g_H–O(r) plots, where the first maximum decreases from a distance of 1.875 Å at the GS to 1.725 Å at the L_a state while staying the same as the GS for the L_b state. The peak varies for GS to L_b and L_a from g_H–O(r) values of 1.259–1.400 and 1.637, and the coordination number is 0.864, 1.028, and 1.051, respectively. Based on the coordination number, only one water molecule HB with indole is in the ground and two excited states (with a difference of 0.15 Å on the HB length between L_b and L_a). Quantum mechanical calculations agree with the classical MD simulations that the HB interaction between indole and water is stronger in the L_a state, although the change is not as large as in MD. The ground and excited S₁(L_a) state minima of the indole +1H₂O complex optimized at the CAM-B3LYP+D3/6-311+G(d) level give a HB distance of 1.88 Å and 1.82 Å, respectively.

In order to examine whether MD simulations can model the state inversion, we assessed how much the water stabilizes indole when it is in the ground, L_b and L_a states, by computing the free energy of solvation, ΔG^sol$v$, using FEP/MD simulations (see Table I). The L_b state is predicted by the MD simulation to be less stabilized by the water solvent than the ground state (GS) with a free energy of solvation of −21.23 kJ/mol ± 0.03 kJ/mol and −23.20 kJ/mol ± 0.04 kJ/mol, respectively. The L_a state is found to be stabilized by almost twice the energy of the L_b state with a free energy of solvation of −42.06 kJ/mol ± 0.04 kJ/mol. Without scaling the attractive C₆ interaction of vdW, the electrostatic interaction between indole and the water solvents contributes −38.80 kJ/mol ± 0.03 kJ/mol to stabilize the L_a state. Scaling the parameters further stabilizes the L_a state from the Coulombic potential alone. From looking at the magnitude of the dipole moments of the electronic GS, L_b state, and L_a state of 2.13 D, 2.23 D, and 5.25 D, it is obvious that the L_b state should be more stabilized by solvation, as is correctly predicted. The difference in the free energy of solvation for the L_b and L_a excited states is a reflection of parameterizing the Coulombic partial charges at the EOM-CCSD/6-311+G(d) level of theory for each state.

The difference between ΔG^sol$v$ for the ground state and the excited state minimum reflects the amount of energy it takes for water to reorganize itself to the electrostatics of indole upon excited state relaxation in aqueous solution, i.e., ΔG^reorg (see Table I). The FEP/MD simulations predict that it takes +1.97 kJ/mol ± 0.05 kJ/mol and −18.86 kJ/mol ± 0.06 kJ/mol for water to respond to the change in the electrostatics of the ground state and reorient themselves to the electrostatics of the excited L_b and L_a states. The positive sign for $ΔGLbreorg$ reflects that the reorganization of waters from the ground state to the L_b minimum is not a favorable process, whereas the negative sign for $ΔGLareorg$ is.

We add the free energy of solvent reorganization to the quantum mechanically computed adiabatic free energy $ΔGgasadiab$ to determine whether the minima of the excited state potential energy surfaces become redshifted or blueshifted in an explicit solvation model using the FEP/MD simulation approach with a non-polarizable force field. The results are shown in Fig. 10. In our previous work, we computed the adiabatic energy from the ground state to the L_b and L_a minima with EOM-CCSD/6-311+G(d) in the gas phase to be 4.70 eV and 5.04 eV, respectively. When adding the zero point energies and thermal effects, the adiabatic free energies are found to be 4.58 eV and 4.86 eV (shown in Fig. 10), respectively. Adding $ΔGLbreorg$ and $ΔGLareorg$ to the gas phase adiabatic free energies, we get that the adiabatic free energies for indole in solution are 4.60 eV and 4.66 eV for L_b and L_a, respectively. Although the L_a state became largely red-shifted, it still lies higher in energy than the L_b state by +0.06 eV. This indicates that the electrostatic effects and explicit solute–solvent interactions are not enough to reproduce this unique energy level inversion at emission. On the other hand, the EOM-CCSD/PCM predicts that L_a becomes barely stabilized past the L_b minimum by a small difference of −0.0018 eV. This indicates that polarization of the solute by the solvent is a crucial factor for observing this special PES inversion of the excited ππ^* states of indole at fluorescence in aqueous solution.

In order to better understand how the different solvation effects (such as polarization effects, HB interaction, and explicit solute–solvent and long-range interaction) contribute to the large redshift resulting in an energy level inversion and an increase in fluorescence intensity, we compare the emission maxima and corresponding oscillator strengths from various models to the experimental ones in Fig. 11.

EOM-CCSD accurately predicts the character of the emitting S₁ state in the gas phase to be L_b and in aqueous solution using the state-specific implicit solvation model, PCM-PTES, to be L_a but overestimates both energies by 0.34 eV and 0.70 eV, respectively. The much larger error for L_a reflects both the fact that the PES of this state is more challenging to describe and the larger errors in PCM solvation. Inclusion of one water molecule hydrogen bonded to N–H of indole modeled quantum mechanically and with EFP already switches the character of the emitting state, predicting it to have L_a character, but overestimates the emission energy by 1.04 eV and 1.06 eV, respectively. Although the MD simulations cannot properly predict the inversion (as shown in Sec. III C 2), we can, nevertheless, use snapshots from the L_a simulations to calculate the emission energies at the L_a minimum. When using snapshots from the MD simulations and water molecules within a solvation shell of 4 Å, the average S₁ energy is 4.47 eV with a standard deviation of 0.17 eV, giving an error of 0.9 eV compared to the experiment. Therefore, the error actually increased compared to the PCM.

These results indicate that indole is a very polarizable molecule, in which polarization from the PCM or HB interaction from just one water molecule at the N–H site of indole is enough to stabilize the L_a excited state to becoming the S₁ minimum. Although this specific interaction of 1H₂O at the quantum mechanical level predicted the emitting S₁ state to be L_a, the stabilization is not at the same magnitude as in the PCM. The PCM predicts the vertical emission energy to be 4.29 eV, whereas the inclusion of 1H₂O gives an energy of 4.64 eV. This suggests that polarization effects at the emission are stronger than those in the HB interaction. Modeling 1H₂O with EFP gives the very same results as the full quantum mechanical description at the EOM-CCSD level, confirming once again that using EOM-CCSD/EFP is an accurate and more computationally affordable technique when used to predict the vertical emission of a solute with a larger solvation shell. The oscillator strength (shown in Fig. 11) is more sensitive to both HB and polarization, and it appears that both are needed to capture the intensity.

It is challenging to model fluorescence in explicit solvation. The models used here to model the fluorescence of indole in aqueous solution are not as extensive as the ones we used for absorption. We used a simple classical non-polarizable force field and molecular dynamics of indole in the L_a excited state, and the bonds were constrained, which means that the vibrational effects are not included. Adding vibrational effects would require reparameterization of all the parameters in the force field, which is a complicated task. In addition, the force field uses fixed charges, therefore not including polarization effects as well (so, equivalent to the C-MM model in absorption). Thus, only explicit solute–solvent electrostatic interactions and solvation dynamical effects are involved. We compare how modeling the electrostatic effects from this explicit solvent model, where HB is prominent, reproduces the fluorescence spectrum of indole in aqueous solution in comparison to that of the implicit PCM. Figure 12 shows the spectra produced from the MD simulation compared to the one produced from the PCM and the experimental spectrum.

Experimental emission λ_max is measured at 3.54 eV. As discussed before, EOM-CCSD has an error when describing the excited states of indole, and since here we are interested in the solvent effect, we want to subtract that error. We will use the same error as in absorption (0.46 eV according to our previous work), although the error may be larger in emission.⁴³ Fluorescence spectra simulated from the classical L_a-MM/MD model predict λ_max to be at 4.03 eV (when corrected), and PCM-PTES predicts it to be at 3.83 eV (when corrected). The two methods differ in the predicted emission maximum by 0.2 eV. Interestingly, the PTES-PCM is closer to the experimental emission maximum compared to the MD results. This shows the importance of mutual polarization, which is present in the PCM.

The modeling for fluorescence in this work is limited. We offer two suggestions for future considerations to indirectly introduce polarization effects in modeling the fluorescence spectra in aqueous solution. One is polarizing the excited state density by computing the fixed Coulombic charges of the excited state using a polarizable continuum model. The second is either scaling the non-bonding parameters of the AMOEBA force field or using the partial charges of the emitting exited state in the polarizable force field to model the solute in the excited state. The second case would be a challenge of its own, considering the sophistication of the polarizable force field in comparison to the classical non-polarizable force field. Despite its shortcomings, the PCM is still a reasonable approach to calculate emission spectra and is the only practical way to calculate stationary points on the excited state potential energy surface of molecules in solution.

In conclusion, we explored the solvation effects on the absorption and fluorescence spectra of indole using various implicit and explicit solvation models. We find that the simulation of a system composed of 394 atoms using ab initio Car–Parrinello MD needs to be longer than 8.3 ps to produce uncorrelated dynamical configurations, which would be computationally exhaustive. Although the spectrum reproduces the position of the peaks, it predicts a much narrower spectrum compared to the experiment. We identify a polarizable AMOEBA force field and molecular dynamics as an affordable theoretical method for capturing indole–water configurations that can reproduce the absorption spectra of indole in solution with some level of accuracy. Modeling the mutual polarization effects between the solute and solvent molecules in the ground state was found necessary to reproduce the different features of the absorption spectrum in solution. SCAN-CPMD produces the best agreement on the position of the peaks, while AMOEBA produces the best line shape. We also show that the lack of polarization effects in the non-polarizable counterparts in comparison to the PCM renders them incapable of predicting the unique energy level inversion between the L_b and L_a excited states that indole exhibits at emission upon solvation. Although using the non-polarizable force field with molecular dynamics did not accurately reproduce the fluorescence spectrum, it portrayed a simple and inexpensive way to overcome many of the challenges that calculating excited state relaxation in explicit solvent would generally present. We are optimistic that this study sheds light into the spectra of indole in aqueous solution and provides an insight into the important factors involved in modeling solvated molecules.

See the supplementary material for additional information on methodology, benchmarking, and results.

S.M. and S.A. acknowledge the support from the National Science Foundation under Grant No. CHE-1800171. S.A. also thanks Temple University for a First Summers Research Initiative Award. This research includes calculations carried out on HPC resources supported, in part, by the National Science Foundation through major research instrumentation (Grant No. 1625061) and by the US Army Research Laboratory under Contract No. W911NF-16-2-0189. V.C. acknowledges the support from the National Institutes of Health through Grant No. GM131048.

The data that support the findings of this study are available from the corresponding author upon reasonable request.