Modeling the optical spectra of molecules in solution presents a challenge, so it is important to understand which of the solvation effects (i.e., electrostatics, mutual polarization, and hydrogen bonding interactions between solute and solvent molecules) are crucial in 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. In this study, we use various implicit and explicit solvation models, such as molecular dynamics coupled with non-polarizable and polarizable force fields, as well as Car–Parrinello molecular dynamics, to model the absorption and fluorescence spectra of indole in aqueous solution. The excited states are computed using the equation of motion coupled cluster with single and double excitations combined with the effective fragment potential to represent water molecules, which we found to be a computationally efficient approach for modeling large solute–solvent clusters at a high level of quantum theory. We find that modeling mutual polarization, compared to other solvation effects, is a dominating factor for accurately reproducing the position of the peaks and spectral line shape of the absorption spectrum of indole in solution. We present an in-depth analysis of the influence that different solvation models have on the electronic excited states responsible for the features of the absorption spectra. Modeling fluorescence is more challenging since it is hard to reproduce even the correct emitting state, and force field parameters need to be re-evaluated.

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.13 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).1 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.1 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.13 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 S1 are calculated in the presence of the solvent.12 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.20 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.12 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.22 

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),23 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.20 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.26 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.19 

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.

1. Polarizable continuum model (PCM)

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−15 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 S1 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−9 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 S1 state of the solute, and the energy is resolved using the equilibrium solvation approach.

2. Indole–water dimer

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

FIG. 1.

Indole +1H2O complex at the ground state minimum.

FIG. 1.

Indole +1H2O complex at the ground state minimum.

Close modal

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 S1 minimum from which fluorescence occurs. All geometry optimization and vibrational frequency calculations were performed on the Gaussian16 software.77 A single point calculation was performed at the optimized S1 minimum of the indole +1H2O 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 +1H2O complex. Indole was described as the QM region computed at the EOM-CCSD level, and the 1H2O molecule was modeled with EFP with parameters computed at the same level of theory.

1. Non-polarizable force field

A structural model of indole solvated by a cubic box of 3291 water molecules was simulated using Gromacs 2016.3.78 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 approach79 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.80 We used the SPC model for water molecules.81 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.82 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 algorithms83,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,85 and the equation of motion was integrated using the leapfrog86 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.

2. AMOEBA polarizable force field

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.87 Parameters from amoeba0988 and amoeba0389 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.90 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.

3. Car–Parrinello

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.92 Norm-conserving HSCV-PBE pseudopotentials93 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.

FIG. 2.

Random snapshot showing indole and water molecules at a distance of (a) 4 Å, (b) 6 Å, and (c) 8 Å from indole in UC-MM/MD.

FIG. 2.

Random snapshot showing indole and water molecules at a distance of (a) 4 Å, (b) 6 Å, and (c) 8 Å from indole in UC-MM/MD.

Close modal

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.94 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 Lb and La 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 Lb and La states of indole, respectively.

The Coulombic potential was parameterized by computing the partial charges of an isolated indole molecule using the densities of La and Lb 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 Lb and La 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 Lb and La. 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,

(1)

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 C12(i, j) and C6(i, j). Here, we consider how the electronic excitation affects these parameters and, in particular, C6(i, j), which according to the Slater–Kirkwood equation101 depends on atomic polarizability as follows:

(2)

where E0 is the ionization energy, a0 is the atomic radius, αi is the polarizability (in units of V−1), NA is the Avogadro number, and Neff,i is the number of effective electrons in atom i. Since the atomic polarizability αi is directly proportional to the atomic volume,21 the denominator in the right-hand side of Eq. (2) is proportional to ViNeff,i, where Vi is the atomic volume, and it is, therefore, also proportional to 1ρi, where ρi is the atomic density,

(3)

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

(4)

by a factor k equal to

(5)

where C6s, Vs, and ρs are the C6 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.102 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 Lb excited state. Upon solvation, the La excited state is redshifted below the Lb excited state, becoming the S1 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 Lb and La 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).

FIG. 3.

Diagram of the potential energy surface of the electronic ground state (gs) and excited state (es) in the gas phase (black solid curve) and in solution (black dotted curve). The diagram portrays the process to compute the adiabatic free energy in an explicit solvent (blue dotted arrow) in three steps. (i) The desolvation of the solute in the ground state (yellow arrow), (ii) electronic excitation (blue solid arrow) followed by excited state relaxation, and (iii) solvation of the solute in the excited state (red arrow).

FIG. 3.

Diagram of the potential energy surface of the electronic ground state (gs) and excited state (es) in the gas phase (black solid curve) and in solution (black dotted curve). The diagram portrays the process to compute the adiabatic free energy in an explicit solvent (blue dotted arrow) in three steps. (i) The desolvation of the solute in the ground state (yellow arrow), (ii) electronic excitation (blue solid arrow) followed by excited state relaxation, and (iii) solvation of the solute in the excited state (red arrow).

Close modal

To evaluate whether modeling the Lb and La excited states of indole in explicit solvation with the Lb-MM/MD model and La-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, ΔGreorg. We define ΔGreorg as follows:

(6)

We compute the free energy of solvation, ΔGsolv, starting with the (GS) C-MM/MD, Lb-MM/MD, and La-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 La-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 La 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 S1 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 gH–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 gH–O(r) is not available, we cannot characterize the accuracy of each method. Instead, we discuss similarities and differences between different theoretical models. gH–O(r) are plotted in Fig. 4.

FIG. 4.

Radial distribution function (RDF) between amino hydrogen (H) of indole and water oxygen (O) obtained through sampling of the ground state ensemble with the classical non-polarizable force field (UC-MM/MD), polarizable force field (AMOEBA MM/MD), and plane wave-pseudopotential SCAN-DFT Car–Parrinello MD. The hydrogen–oxygen distance is scanned in increments of 0.15 Å. The RDF obtained from the BYLP-BOMD simulation is from Ref. 103.

FIG. 4.

Radial distribution function (RDF) between amino hydrogen (H) of indole and water oxygen (O) obtained through sampling of the ground state ensemble with the classical non-polarizable force field (UC-MM/MD), polarizable force field (AMOEBA MM/MD), and plane wave-pseudopotential SCAN-DFT Car–Parrinello MD. The hydrogen–oxygen distance is scanned in increments of 0.15 Å. The RDF obtained from the BYLP-BOMD simulation is from Ref. 103.

Close modal

1. First solvation shell

The peak in gH–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 +1H2O 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 +1H2O complex and the peak of the gH–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 AMOEBA90 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, gH–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. gH–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 +1H2O complex with that predicted from the MD simulations.

2. Second solvation shell

The UC and AMOEBA MM/MD simulations do not show pronounced second peaks in gH–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 gH–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. gH–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 gO–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 gH–O(r), and as a result, some of features of gH–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.

FIG. 5.

Absorption spectra simulated from four excited states of 50 uncorrelated indole–water configurations sampled with various solvation models and calculated at the EOM-CCSD/EFP level of theory. A Lorentzian line shape function with a linewidth of 0.2 eV was used to fit the electronic energy distribution. Simulated spectra are overlaid with the experimental spectra from Ref. 104. The intensities of the peaks have been normalized to that of observed λmax. In (a), the theoretical spectra are shifted by 0.46 eV to shorter energies to correct for the gas phase EOM-CCSD error. In (b), they are all shifted so that their maxima coincide.

FIG. 5.

Absorption spectra simulated from four excited states of 50 uncorrelated indole–water configurations sampled with various solvation models and calculated at the EOM-CCSD/EFP level of theory. A Lorentzian line shape function with a linewidth of 0.2 eV was used to fit the electronic energy distribution. Simulated spectra are overlaid with the experimental spectra from Ref. 104. The intensities of the peaks have been normalized to that of observed λmax. In (a), the theoretical spectra are shifted by 0.46 eV to shorter energies to correct for the gas phase EOM-CCSD error. In (b), they are all shifted so that their maxima coincide.

Close modal

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.

1. λMAX peak

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.43 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.

2. Second peak

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.

3. Spectral shape

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.

4. A closer look

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 La and Lb. 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 Lb, La, 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.

FIG. 6.

Absorption spectra calculated separately for the ππ* Lb and La excited states and the two πσ* states. Top row [panels (a)–(c)]: Excited states calculated using indole and water molecules are used. Bottom row [panel (d)–(f)]: Excited states were calculated using only indole in order to separate the vibrational effects from solvation interactions.

FIG. 6.

Absorption spectra calculated separately for the ππ* Lb and La excited states and the two πσ* states. Top row [panels (a)–(c)]: Excited states calculated using indole and water molecules are used. Bottom row [panel (d)–(f)]: Excited states were calculated using only indole in order to separate the vibrational effects from solvation interactions.

Close modal
a. Electronic states.

It is obvious from Fig. 6 that the two main peaks seen in all spectra arise from the electronic states La and Lb. The higher energy peak that corresponds to λmax is due to the La state, while the second peak at lower energy is due to the Lb state (mostly). It is also clear that while the Lb band has a single maximum and is described well by a Lorentzian, this is not the case for the La 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 La and the Lb states since one of the subbands in La blends with the Lb band. The combined band is broader, and the maximum is shifted to higher energies compared to the Lb 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 La, 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.

FIG. 7.

Absorption spectra simulated with (solid line) and without (dashed line) the inclusion of πσ* states.

FIG. 7.

Absorption spectra simulated with (solid line) and without (dashed line) the inclusion of πσ* states.

Close modal
b. Vibrational effects on the electronic bands.

In order to explore the origin of the detailed structure in the La 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 La 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−1 from the peak [panel (a) of Fig. 6]. In the La spectra simulated from the AMOEBA MM/MD configurations [panel (b) of Fig. 6], the separation is 1033 cm−1. SCAN-CPMD, on the other hand, only predicts a single peak [panel (c) of Fig. 6]. It is clear from these spectra that the La band exhibits a vibrational progression, while Lb does not. This is in agreement with their minimum energy geometries and how they differ from the ground state minimum. The Lb minimum is very similar to the ground state minimum, so the Franck–Condon factors do not favor a vibrational progression. The La 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.

c. Solvation effects on the electronic bands.

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 La and Lb peaks. It also increases the separation of the vibrational progressions in La. In UC-MM, the separation between the two peaks is 1722 cm−1, while in AMOEBA, it is 1301 cm−1. In AMOEBA, the gap between the La and Lb peaks decreases with solvation to the extent that Lb becomes almost completely hidden beneath the La 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.

5. Absorption spectra with TDDFT

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 Lb(1ππ*) state to be brighter than La(2ππ*). Modeling the solvent explicitly, on the other hand, predicts that indole absorbs mostly to the S1 excited state and not S2. It also predicts the wrong ordering of states at the ground state minimum, with the S1 state being La and the S2 excited state being Lb. 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.

FIG. 8.

Absorption spectra simulated from four excited states and oscillator strengths of 50 configurations of indole clusters sampled with UC-MM/MD, AMOEBA MM/MD, and SCAN-CPMD simulations, computed at the CAM-B3LYP+D3/6-311+G(d) level of theory (tabulated in Tables S9 and S10) using a linewidth of 0.2 eV.

FIG. 8.

Absorption spectra simulated from four excited states and oscillator strengths of 50 configurations of indole clusters sampled with UC-MM/MD, AMOEBA MM/MD, and SCAN-CPMD simulations, computed at the CAM-B3LYP+D3/6-311+G(d) level of theory (tabulated in Tables S9 and S10) using a linewidth of 0.2 eV.

Close modal

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.39 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 S1 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 La state is found to be largely stabilized by the water solvent to the extent that an energy inversion between the La and Lb excited state takes place at the fluorescence where La becomes the emitting state.40,48,104–106 Hence, experimentally, it is observed that while in the gas phase the emission occurs from the Lb state, the ordering of the states is inverted in aqueous solution and emission occurs from La. 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.

1. Excited state sampling

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 Lb and La 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 Lb-MM/MD and La-MM/MD. The changes in the solvation structure upon excitation are shown by plotting gH–O(r) for the ground and excited Lb and La states (shown in Fig. 9).

FIG. 9.

Radial distribution function of H bonded to N of indole and gH–O(r) in the excited Lb and La states to O of water modeled with C-MM/MD [bin 0.15].

FIG. 9.

Radial distribution function of H bonded to N of indole and gH–O(r) in the excited Lb and La states to O of water modeled with C-MM/MD [bin 0.15].

Close modal

The HB strength with water becomes stronger in the La state, as is indicated in the gH–O(r) plots, where the first maximum decreases from a distance of 1.875 Å at the GS to 1.725 Å at the La state while staying the same as the GS for the Lb state. The peak varies for GS to Lb and La from gH–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 Lb and La). Quantum mechanical calculations agree with the classical MD simulations that the HB interaction between indole and water is stronger in the La state, although the change is not as large as in MD. The ground and excited S1(La) state minima of the indole +1H2O complex optimized at the CAM-B3LYP+D3/6-311+G(d) level give a HB distance of 1.88 Å and 1.82 Å, respectively.

2. Examining energy level inversion using FEP/MD

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, Lb and La states, by computing the free energy of solvation, ΔGsolv, using FEP/MD simulations (see Table I). The Lb 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 La state is found to be stabilized by almost twice the energy of the Lb state with a free energy of solvation of −42.06 kJ/mol ± 0.04 kJ/mol. Without scaling the attractive C6 interaction of vdW, the electrostatic interaction between indole and the water solvents contributes −38.80 kJ/mol ± 0.03 kJ/mol to stabilize the La state. Scaling the parameters further stabilizes the La state from the Coulombic potential alone. From looking at the magnitude of the dipole moments of the electronic GS, Lb state, and La state of 2.13 D, 2.23 D, and 5.25 D, it is obvious that the Lb state should be more stabilized by solvation, as is correctly predicted. The difference in the free energy of solvation for the Lb and La 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.

TABLE I.

Free energy of solvation, ΔGsolv, and reorganization energy, ΔGreorg, for GS and excited states, Lb and La.

C-MM/MD modelΔGsolv (kJ/mol)ΔGreorg (kJ/mol)
GS −23.20 ± 0.04  
Lb −21.23 ± 0.03 +1.97 ± 0.05 
La −42.06 ± 0.04 −18.86 ± 0.06 
C-MM/MD modelΔGsolv (kJ/mol)ΔGreorg (kJ/mol)
GS −23.20 ± 0.04  
Lb −21.23 ± 0.03 +1.97 ± 0.05 
La −42.06 ± 0.04 −18.86 ± 0.06 

The difference between ΔGsolv 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., ΔGreorg (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 Lb and La states. The positive sign for ΔGLbreorg reflects that the reorganization of waters from the ground state to the Lb 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 Lb and La 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 Lb and La, respectively. Although the La state became largely red-shifted, it still lies higher in energy than the Lb 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 La becomes barely stabilized past the Lb 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.

FIG. 10.

Potential energy surface diagram of the ground (gs) and the excited states, Lb and La, portraying the effect of solvation (dotted line) on indole in the gas phase surface (solid line). ΔGsolv for each state is indicated, as well as the adiabatic energies ΔGadiab in the gas phase (ΔGgasadiab) and solution (ΔGsolnadiab).

FIG. 10.

Potential energy surface diagram of the ground (gs) and the excited states, Lb and La, portraying the effect of solvation (dotted line) on indole in the gas phase surface (solid line). ΔGsolv for each state is indicated, as well as the adiabatic energies ΔGadiab in the gas phase (ΔGgasadiab) and solution (ΔGsolnadiab).

Close modal

3. The S1 minimum

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.

FIG. 11.

(a) Vertical emission energies and (b) oscillator strengths computed using different solvation models and their corresponding fluorescence intensity (oscillator strength).

FIG. 11.

(a) Vertical emission energies and (b) oscillator strengths computed using different solvation models and their corresponding fluorescence intensity (oscillator strength).

Close modal

EOM-CCSD accurately predicts the character of the emitting S1 state in the gas phase to be Lb and in aqueous solution using the state-specific implicit solvation model, PCM-PTES, to be La but overestimates both energies by 0.34 eV and 0.70 eV, respectively. The much larger error for La 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 La 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 La simulations to calculate the emission energies at the La minimum. When using snapshots from the MD simulations and water molecules within a solvation shell of 4 Å, the average S1 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 La excited state to becoming the S1 minimum. Although this specific interaction of 1H2O at the quantum mechanical level predicted the emitting S1 state to be La, 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 1H2O gives an energy of 4.64 eV. This suggests that polarization effects at the emission are stronger than those in the HB interaction. Modeling 1H2O 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.

4. Fluorescence spectra

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 La 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.

FIG. 12.

Simulated fluorescence spectra from the C-MM/MD model of indole in the La excited state and PCM-PTES shifted by 0.46 eV and overlaid with the experimental emission spectra taken from Ref. 104.

FIG. 12.

Simulated fluorescence spectra from the C-MM/MD model of indole in the La excited state and PCM-PTES shifted by 0.46 eV and overlaid with the experimental emission spectra taken from Ref. 104.

Close modal

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.43 Fluorescence spectra simulated from the classical La-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 Lb and La 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.

1.
T.
Giovannini
,
F.
Egidi
, and
C.
Cappelli
, “
Molecular spectroscopy of aqueous solutions: A theoretical perspective
,”
Chem. Soc. Rev.
49
,
5664
5677
(
2020
).
2.
M. S.
de Vries
and
P.
Hobza
, “
Gas-phase spectroscopy of biomolecular building blocks
,”
Annu. Rev. Phys. Chem.
58
,
585
612
(
2007
).
3.
M.
Schreiber
,
M. R.
Silva
, Jr.
,
S. P. A.
Sauer
, and
W.
Thiel
, “
Benchmarks for electronically excited states: CASPT2, CC2, CCSD, and CC3
,”
J. Chem. Phys.
128
,
134110
(
2008
).
4.
R.
Improta
,
F.
Santoro
, and
L.
Blancafort
, “
Quantum mechanical studies on the photophysics and the photochemistry of nucleic acids and nucleobases
,”
Chem. Rev.
116
,
3540
3593
(
2016
).
5.
L.
González
,
D.
Escudero
, and
L.
Serrano-Andrés
, “
Progress and challenges in the calculation of electronic excited states
,”
ChemPhysChem
13
,
28
51
(
2012
).
6.
G.-Y.
Li
and
K.-L.
Han
, “
The sensing mechanism studies of the fluorescent probes with electronically excited state calculations
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1351
(
2018
).
7.
J.
Linnanto
and
J.
Korppi-Tommola
, “
Quantum chemical simulation of excited states of chlorophylls, bacteriochlorophylls and their complexes
,”
Phys. Chem. Chem. Phys.
8
,
663
687
(
2006
).
8.
M.
Orozco
and
F. J.
Luque
, “
Theoretical methods for the description of the solvent effect in biomolecular systems
,”
Chem. Rev.
100
,
4187
4225
(
2000
).
9.
J.
Tomasi
and
M.
Persico
, “
Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent
,”
Chem. Rev.
94
,
2027
(
1994
).
10.
J.
Tomasi
,
B.
Mennucci
, and
R.
Cammi
, “
Quantum mechanical continuum solvation models
,”
Chem. Rev.
105
,
2999
(
2005
).
11.
R.
Improta
and
V.
Barone
, “
Excited states behavior of nucleobases in solution: Insights from computational studies
,” in
Photoinduced Phenomena in Nucleic Acids I: Nucleobases in the Gas Phase and in Solvents
, Topics in Current Chemistry-Series, edited by
M.
Barbatti
,
A. C.
Borin
, and
S.
Ullrich
(
Springer
,
2015
), Vol. 355, pp.
329
357
.
12.
T. J.
Zuehlsdorff
and
C. M.
Isborn
, “
Modeling absorption spectra of molecules in solution
,”
Int. J. Quantum Chem.
119
,
e25719
(
2019
).
13.
N. H.
List
,
J. M. H.
Olsen
, and
J.
Kongsted
, “
Excited states in large molecular systems through polarizable embedding
,”
Phys. Chem. Chem. Phys.
18
,
20234
20250
(
2016
).
14.
D.
Loco
,
N.
Gelfand
,
S.
Jurinovich
,
S.
Protti
,
A.
Mezzetti
, and
B.
Mennucci
, “
Polarizable QM/classical approaches for the modeling of solvation effects on UV–vis and fluorescence spectra: An integrated strategy
,”
J. Phys. Chem. A
122
,
390
397
(
2018
).
15.
T.
Petrenko
and
F.
Neese
, “
Analysis and prediction of absorption band shapes, fluorescence band shapes, resonance Raman intensities, and excitation profiles using the time-dependent theory of electronic spectroscopy
,”
J. Chem. Phys.
127
,
164319
(
2007
).
16.
B.
Mennucci
, “
Polarizable continuum model
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
386
404
(
2012
).
17.
D.
Loco
and
L.
Cupellini
, “
Modeling the absorption lineshape of embedded systems from molecular dynamics: A tutorial review
,”
Int. J. Quantum Chem.
119
,
e25726
(
2019
).
18.
J. J.
Nogueira
and
L.
González
, “
Computational photophysics in the presence of an environment
,”
Annu. Rev. Phys. Chem.
69
,
473
497
(
2018
).
19.
A.
DeFusco
,
N.
Minezawa
,
L. V.
Slipchenko
,
F.
Zahariev
, and
M. S.
Gordon
, “
Modeling solvent effects on electronic excited states
,”
J. Phys. Chem. Lett.
2
,
2184
2192
(
2011
).
20.
S.
Gómez
,
T.
Giovannini
, and
C.
Cappelli
, “
Absorption spectra of xanthines in aqueous solution: A computational study
,”
Phys. Chem. Chem. Phys.
22
,
5929
5941
(
2020
).
21.
P.
Politzer
,
P.
Jin
, and
J. S.
Murray
, “
Atomic polarizability, volume and ionization energy
,”
J. Chem. Phys.
117
,
8197
(
2002
).
22.
V.
Barone
,
J.
Bloino
,
M.
Biczysko
, and
F.
Santoro
, “
Fully integrated approach to compute vibrationally resolved optical spectra: From small molecules to macrosystems
,”
J. Chem. Theory Comput.
5
,
540
554
(
2009
).
23.
P.
Cieplak
,
F.-Y.
Dupradeau
,
Y.
Duan
, and
J.
Wang
, “
Polarization effects in molecular mechanical force fields
,”
J. Phys.: Condens. Matter
21
,
333102
(
2009
).
24.
S.
Patel
and
C. L.
Brooks
 III
, “
Fluctuating charge force fields: Recent developments and applications from small molecules to macromolecular biological systems
,”
Mol. Simul.
32
,
231
249
(
2006
).
25.
T.
Giovannini
,
M.
Macchiagodena
,
M.
Ambrosetti
,
A.
Puglisi
,
P.
Lafiosca
,
G.
Lo Gerfo
,
F.
Egidi
, and
C.
Cappelli
, “
Simulating vertical excitation energies of solvated dyes: From continuum to polarizable discrete modeling
,”
Int. J. Quantum Chem.
119
,
e25684
(
2019
).
26.
A.
Puglisi
,
T.
Giovannini
,
L.
Antonov
, and
C.
Cappelli
, “
Interplay between conformational and solvent effects in UV-visible absorption spectra: Curcumin tautomers as a case study
,”
Phys. Chem. Chem. Phys.
21
,
15504
15514
(
2019
).
27.
M. A.
Thompson
, “
QM/MMPOL: A consistent model for solute/solvent polarization. Application to the aqueous solvation and spectroscopy of formaldehyde, acetaldehyde, and acetone
,”
J. Phys. Chem.
100
,
14492
14507
(
1996
).
28.
J. M.
Olsen
,
K.
Aidas
, and
J.
Kongsted
, “
Excited states in solution through polarizable embedding
,”
J. Chem. Theory Comput.
6
,
3721
3734
(
2010
).
29.
R.
Guareschi
,
H.
Zulfikri
,
C.
Daday
,
F. M.
Floris
,
C.
Amovilli
,
B.
Mennucci
, and
C.
Filippi
, “
Introducing QMC/MMpol: Quantum Monte Carlo in polarizable force fields for excited states
,”
J. Chem. Theory Comput.
12
,
1674
1683
(
2016
).
30.
D.
Loco
,
É.
Polack
,
S.
Caprasecca
,
L.
Lagardère
,
F.
Lipparini
,
J.-P.
Piquemal
, and
B.
Mennucci
, “
A QM/MM approach using the AMOEBA polarizable embedding: From ground state energies to electronic excitations
,”
J. Chem. Theory Comput.
12
,
3654
3661
(
2016
).
31.
T.
Giovannini
,
A.
Puglisi
,
M.
Ambrosetti
, and
C.
Cappelli
, “
Polarizable QM/MM approach with fluctuating charges and fluctuating dipoles: The QM/FQFμ model
,”
J. Chem. Theory Comput.
15
,
2233
2245
(
2019
).
32.
M. S.
Gordon
,
M. A.
Freitag
, and
P.
Bandyopadhyay
, “
The effective fragment potential method: A QM-based MM approach to modeling environmental effects in Chemistry
,”
J. Phys. Chem. A
105
,
293
307
(
2001
).
33.
M. S.
Gordon
,
L.
Slipchenko
,
H.
Li
, and
J. H.
Jensen
, “
The effective fragment potential: A general method for predicting intermolecular interactions
,”
Annu. Rep. Comput. Chem.
3
,
177
(
2007
).
34.
M. S.
Gordon
,
D. G.
Fedorov
,
S. R.
Pruitt
, and
L. V.
Slipchenko
, “
Fragmentation methods: A route to accurate calculations on large systems
,”
Chem. Rev.
112
,
632
(
2012
).
35.
M. S.
Gordon
,
Q. A.
Smith
,
P.
Xu
, and
L. V.
Slipchenko
, “
Accurate first principles model potentials for intermolecular interactions
,”
Ann. Rev. Phys. Chem.
64
,
553
(
2013
).
36.
R. J.
Melander
,
M. J.
Minvielle
, and
C.
Melander
, “
Controlling bacterial behavior with indole-containing natural products and derivatives
,”
Tetrahedron
70
,
6363
6372
(
2014
).
37.
J.-H.
Lee
and
J.
Lee
, “
Indole as an intercellular signal in microbial communities
,”
FEMS Microbiol. Rev.
34
,
426
444
(
2010
).
38.
R.
Benjamins
and
B.
Scheres
, “
Auxin: The looping star in plant development
,”
Annu. Rev. Plant Biol.
59
,
443
465
(
2008
).
39.
S.
Arulmozhiraja
and
M. L.
Coote
, “
1La and 1Lb states of indole and azaindole: Is density functional theory inadequate?
,”
J. Chem. Theory Comput.
8
,
575
584
(
2012
).
40.
D.
Brisker-klaiman
and
A.
Dreuw
, “
Explaining level inversion of the La and Lb states of indole and indole derivatives in polar solvents
,”
ChemPhysChem
16
,
1695
1702
(
2015
).
41.
L. S.
Slater
and
P. R.
Callis
, “
Molecular orbital theory of the 1Lb and 1La states of indole
,”
J. Phys. Chem.
99
,
8572
8581
(
1995
).
42.
L.
Serrano-Andrés
and
B. O.
Roos
, “
Theoretical study of the absorption and emission spectra of indole in the gas phase and in a solvent
,”
J. Am. Chem. Soc.
118
,
185
195
(
1996
).
43.
S.
Abou-Hatab
and
S.
Matsika
, “
Theoretical investigation of positional substitution and solvent on n-cyanoindole fluorescent probes
,”
J. Phys. Chem. B
123
,
7424
7435
(
2019
).
44.
N.
Sharma
,
S. K.
Jain
, and
R. C.
Rastogi
, “
Solvatochromic study of excited state dipole moments of some biologically active indoles and tryptamines
,”
Spectrochim. Acta, Part A
66
,
171
176
(
2007
).
45.
M. R.
Eftink
,
L. A.
Selvidge
,
P. R.
Callis
, and
A. A.
Rehms
, “
Photophysics of indole derivatives: Experimental resolution of La and lb transitions and comparison with theory
,”
J. Phys. Chem.
94
,
3469
3479
(
1990
).
46.
B.
Albinsson
and
B.
Nordén
, “
Excited-state properties of the indole chromophore. Electronic transition moment directions from linear dichroism measurements: Effect of methyl and methoxy substituents
,”
J. Phys. Chem.
96
,
6204
6212
(
1992
).
47.
P. L.
Muiño
and
P. R.
Callis
, “
Hybrid simulations of solvation effects on electronic spectra: Indoles in water
,”
J. Chem. Phys.
100
,
4093
(
1994
).
48.
X.
Meng
,
T.
Harricharran
, and
L. J.
Juszczak
, “
A spectroscopic survey of substituted indoles reveals consequences of a stabilized lb transition
,”
Photochem. Photobiol.
89
,
40
50
(
2013
).
49.
J.
Wilke
,
M.
Wilke
,
C.
Brand
,
J. D.
Spiegel
,
C. M.
Marian
, and
M.
Schmitt
, “
Modulation of the La/Lb mixing in an indole derivative: A position-dependent study using 4-, 5-, and 6-fluoroindole
,”
J. Phys. Chem. A
121
,
1597
1606
(
2017
).
50.
J.
Wilke
,
M.
Wilke
,
C.
Brand
,
W. L.
Meerts
, and
M.
Schmitt
, “
On the additivity of molecular fragment dipole moments of 5-substituted indole derivatives
,”
ChemPhysChem
17
,
2736
2743
(
2016
).
51.
M. S.
Walker
,
T. W.
Bednar
, and
R.
Lumry
, “
Exciplex formation in the excited state of indole
,”
J. Chem. Phys.
45
,
3455
(
1966
).
52.
M.
Sun
and
P.-S.
Song
, “
Solvent effects on the fluorescent states of indole derivatives-dipole moments
,”
Photochem. Photobiol.
25
,
3
9
(
1977
).
53.
T. J.
Godfrey
,
H.
Yu
,
S.
Ullrich
,
T. J.
Godfrey
,
H.
Yu
, and
S.
Ullrich
, “
Investigation of electronically excited indole relaxation dynamics via photoionization and fragmentation pump-probe spectroscopy
,”
J. Chem. Phys.
141
,
044314
(
2014
).
54.
H. T.
Yu
,
W. J.
Colucci
,
M. L.
Mclaughlin
, and
M. D.
Barkley
, “
Fluorescence quenching in indoles by excited-state proton transfer
,”
J. Am. Chem. Soc.
114
,
8449
8454
(
1992
).
55.
H.
Lami
and
N.
Glasser
, “
Indole’s solvatochromism revisited
,”
J. Chem. Phys.
84
,
597
604
(
1986
).
56.
J.
Cataln
and
J. P.
Cataln
, “
Questioning the photophysical model for the indole chromophore in the light of evidence obtained by controlling the non-specific effect of the medium with 1-chlorobutane as solvent
,”
Phys. Chem. Chem. Phys.
13
,
15022
15030
(
2011
).
57.
J.
Cataln
, “
The first UV absorption band for indole is not due to two simultaneous orthogonal electronic transitions differing in dipole moment
,”
Phys. Chem. Chem. Phys.
17
,
12515
12520
(
2015
).
58.
R. M.
Levy
,
J. D.
Westbrook
,
D. B.
Kitchen
, and
K.
Krogh-jespersen
, “
Solvent effects on the adiabatic free energy difference between the ground and excited states of methylindole in water
,”
J. Phys. Chem.
95
,
6756
6758
(
1991
).
59.
C.
Kang
,
T. M.
Korter
, and
D. W.
Pratt
, “
Experimental measurement of the induced dipole moment of an isolated molecule in its ground and electronically excited states: Indole and indole-H2O
,”
J. Chem. Phys.
122
,
174301
(
2005
).
60.
G.
Kumar
,
A.
Roy
,
R. S.
Mcmullen
,
S.
Kutagulla
, and
S. E.
Bradforth
, “
The influence of aqueous solvent on the electronic structure and non-adiabatic dynamics of indole explored by liquid-jet photoelectron spectroscopy
,”
Faraday Discuss.
212
,
359
381
(
2018
).
61.
A.
Nenov
,
I.
Rivalta
,
S.
Mukamel
, and
M.
Garavelli
, “
Bidimensional electronic spectroscopy on indole in gas phase and in water from first principles
,”
Comput. Theor. Chem.
1040-1041
,
295
303
(
2014
).
62.
G.
Alagona
,
C.
Ghio
, and
S.
Monti
, “
The effect of small substituents on the properties of indole. An ab initio 6-31g* study
,”
J. Mol. Struct.: THEOCHEM
433
,
203
216
(
1998
).
63.
A. L.
Sobolewski
and
W.
Domcke
, “
Photoinduced charge separation in indole–water clusters
,”
Chem. Phys. Lett.
329
,
130
137
(
2000
).
64.
A. L.
Sobolewski
and
W.
Domcke
, “
Computational studies of the photophysics of hydrogen-bonded molecular system
,”
J. Phys. Chem. A
111
,
11725
11735
(
2007
).
65.
J. J.
Aaron
,
A.
Tine
,
C.
Villiers
,
C.
Parkanyi
, and
D.
Bouin
, “
Electronic absorption and fluorescence spectra of indole derivatives. Quantitative treatment of the substituent effects and a theoretical study
,”
Croat. Chem. Acta
56
,
157
168
(
1983
); available at https://hrcak.srce.hr/194220.
66.
J. R.
Lombardi
, “
Solvatochromic shifts: A reconsideration
,”
J. Phys. Chem. A
102
,
2817
3282
(
1998
).
67.
J. R.
Lombardi
, “
Solvatochromic shifts reconsidered: Field-induced mixing in the nonlinear region and application to indole
,”
J. Phys. Chem. A
103
,
6335
6338
(
1999
).
68.
S. R.
Meech
,
D.
Phillips
, and
A. G.
Lee
, “
On the nature of the fluorescent state of methylated indole derivatives
,”
Chem. Phys.
80
,
317
328
(
1983
).
69.
V.
Barone
,
J.
Bloino
,
S.
Monti
,
A.
Pedone
, and
G.
Prampolini
, “
Fluorescence spectra of organic dyes in solution: A time dependent multilevel approach
,”
Phys. Chem. Chem. Phys.
13
,
2160
2166
(
2011
).
70.
M.
Sulpizi
,
U. F.
Röhrig
,
J.
Hutter
, and
U.
Rothlisberger
, “
Optical properties of molecules in solution via hybrid TDDFT/MM simulations
,”
Int. J. Quantum Chem.
101
,
671
682
(
2005
).
71.
M.
D’Abramo1
,
M.
Aschi
, and
A.
Amadei
, “
Theoretical modeling of UV-vis absorption and emission spectra in liquid state systems including vibrational and conformational effects: Explicit treatment of the vibronic transitions
,”
J. Chem. Phys.
140
,
164104
(
2014
).
72.
M.
Aschi
,
V.
Barone
,
B.
Carlotti
,
I.
Daidone
,
F.
Elisei
, and
A.
Amadei
, “
Photoexcitation and relaxation kinetics of molecular systems in solution: Towards a complete in silico model
,”
Phys. Chem. Chem. Phys.
18
,
28919
28931
(
2016
).
73.
C.
García-Iriepa
,
P.
Gosset
,
R.
Berraud-Pache
,
M.
Zemmouche
,
G.
Taupier
,
K. D.
Dorkenoo
,
P.
Didier
,
J.
Léonard
,
N.
Ferré
, and
I.
Navizet
, “
Simulation and analysis of the spectroscopic properties of oxyluciferin and its analogues in water
,”
J. Chem. Theory Comput.
14
,
2117
2126
(
2018
).
74.
M.
Caricato
, “
Absorption and emission spectra of solvated molecules with the EOM–CCSD–PCM method
,”
J. Chem. Theory Comput.
8
,
4494
4502
(
2012
).
75.
M.
Caricato
, “
Exploring potential energy surfaces of electronic excited states in solution with the EOM-CCSD-PCM method
,”
J. Chem. Theory Comput.
8
,
5081
5091
(
2012
).
76.
M.
Caricato
, “
CCSD-PCM: Improving upon the reference reaction field approximation at no cost
,”
J. Chem. Phys.
135
,
074113
(
2011
).
77.
M. J.
Frisch
,
G. W.
Trucks
,
H. B.
Schlegel
,
G. E.
Scuseria
,
M. A.
Robb
,
J. R.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
G. A.
Petersson
,
H.
Nakatsuji
,
X.
Li
,
M.
Caricato
,
A. V.
Marenich
,
J.
Bloino
,
B. G.
Janesko
,
R.
Gomperts
,
B.
Mennucci
,
H. P.
Hratchian
,
J. V.
Ortiz
,
A. F.
Izmaylov
,
J. L.
Sonnenberg
,
D.
Williams-Young
,
F.
Ding
,
F.
Lipparini
,
F.
Egidi
,
J.
Goings
,
B.
Peng
,
A.
Petrone
,
T.
Henderson
,
D.
Ranasinghe
,
V. G.
Zakrzewski
,
J.
Gao
,
N.
Rega
,
G.
Zheng
,
W.
Liang
,
M.
Hada
,
M.
Ehara
,
K.
Toyota
,
R.
Fukuda
,
J.
Hasegawa
,
M.
Ishida
,
T.
Nakajima
,
Y.
Honda
,
O.
Kitao
,
H.
Nakai
,
T.
Vreven
,
K.
Throssell
,
J. A.
Montgomery
, Jr.
,
J. E.
Peralta
,
F.
Ogliaro
,
M. J.
Bearpark
,
J. J.
Heyd
,
E. N.
Brothers
,
K. N.
Kudin
,
V. N.
Staroverov
,
T. A.
Keith
,
R.
Kobayashi
,
J.
Normand
,
K.
Raghavachari
,
A. P.
Rendell
,
J. C.
Burant
,
S. S.
Iyengar
,
J.
Tomasi
,
M.
Cossi
,
J. M.
Millam
,
M.
Klene
,
C.
Adamo
,
R.
Cammi
,
J. W.
Ochterski
,
R. L.
Martin
,
K.
Morokuma
,
O.
Farkas
,
J. B.
Foresman
, and
D. J.
Fox
, Gaussian16 Revision B.01,
Gaussian, Inc.
,
Wallingford, CT,
2016
.
78.
M.
Abraham
,
D.
van der Spoel
,
E.
Lindahl
,
B.
Hess
, and
the GROMACS development team
, Gromacs user manual version 2016.3,
2017
.
79.
C. M.
Breneman
and
K. B.
Wiberg
, “
Determining atom-centered monopoles from molecular electrostatic potentials - the need for high sampling density in formamide conformational-analysis
,”
J. Comput. Chem.
11
,
361
373
(
1990
).
80.
D.
Case
,
D.
Cerutti
,
T.
Cheatham
 III
,
T.
Darden
,
R. E.
Duke
,
T. J.
Giese
,
H.
Gohlke
,
A. W.
Goetz
,
D.
Greene
,
N.
Homeyer
,
S.
Izadi
,
A.
Kovalenko
,
T.
Lee
,
S.
LeGrand
,
P.
Li
,
C.
Lin
,
J.
Liu
,
T.
Luchko
,
R.
Luo
,
D.
Mermelstein
,
K. M.
Merz
,
G.
Monard
,
H.
Nguyen
,
I.
Omelyan
,
A.
Onufriev
,
F.
Pan
,
R.
Qi
,
D.
Roe
,
A.
Roitber
,
C.
Sagui
,
C. L.
Simmerling
,
W. M.
Botello-Smith
,
J.
Swails
,
R.
Walker
,
J.
Wang
,
R. M.
Wolf
,
X.
Wu
,
L.
Xiao
,
D. M.
York
, and
P. A.
Kollman
, Amber 2017,
University of California
,
San Francisco
,
2017
.
81.
H. J. C.
Berendsen
,
J. P. M.
Postma
,
W. F.
van Gunsteren
, and
J.
Hermans
, “
Interaction models for water in relation to protein hydration
,” in , edited by B. Pullman (
Springer, Dordrecht
,
1981
), Vol. 14, pp.
331
342
.
82.
G.
Bussi
,
D.
Donadio
, and
M.
Parrinello
, “
Canonical sampling through velocity rescaling
,”
J. Chem. Phys.
126
,
014101
(
2007
).
83.
M.
Parrinello
and
A.
Rahman
, “
Crystal structure and pair potentials: A molecular-dynamics study
,”
Phys. Rev. Lett.
45
,
1196
1199
(
1980
).
84.
M.
Parrinello
and
A.
Rahman
, “
Polymorphic transitions in single crystals: A new molecular dynamics method
,”
J. Appl. Phys.
52
,
7182
7190
(
1981
).
85.
B.
Hess
,
H.
Bekker
,
H. J. C.
Berendsen
, and
J. G. E. M.
Fraaije
, “
Lincs: A linear constraint solver for molecular simulations
,”
J. Comput. Chem.
18
,
1463
1472
(
1997
).
86.
W. F. V.
Gunsteren
and
H. J. C.
Berendsen
, “
A leap-frog algorithm for stochastic dynamics
,”
Mol. Simul.
1
,
173
185
(
1988
).
87.
J. A.
Rackers
,
Z.
Wang
,
C.
Lu
,
M. L.
Laury
,
L.
Lagardère
,
M. J.
Schnieders
,
J.-P.
Piquemal
,
P.
Ren
, and
J. W.
Ponder
, “
Tinker 8: Software tools for molecular design
,”
J. Chem. Theory Comput.
14
,
5273
5289
(
2018
).
88.
P.
Ren
and
J. W.
Ponder
, “
Consistent treatment of inter- and intramolecular polarization in molecular Mechanics calculations
,”
J. Comput. Chem.
23
,
1497
1506
(
2002
).
89.
P.
Ren
and
J. W.
Ponder
, “
Polarizable atomic multipole water model for molecular Mechanics simulation
,”
J. Phys. Chem. B
107
,
5933
5947
(
2003
).
90.
J. W.
Ponder
,
C.
Wu
,
P.
Ren
,
V. S.
Pande
,
J. D.
Chodera
,
M. J.
Schnieders
,
I.
Haque
,
D. L.
Mobley
,
D. S.
Lambrecht
,
R. A.
Distasio
, Jr.
,
M.
Head-Gordon
,
G. N. I.
Clark
,
M. E.
Johnson
, and
T.
Head-Gordon
, “
Current status of the AMOEBA polarizable force field
,”
J. Phys. Chem. B
114
,
2549
2564
(
2010
).
91.
C.
Pan
,
C.
Liu
,
J.
Peng
,
P.
Ren
, and
X.
Huang
, “
Three-site and five-site fixed-charge water models compatible with AMOEBA force field
,”
J. Comput. Chem.
41
,
1034
1044
(
2020
).
92.
P.
Giannozzi
,
S.
Baroni
,
N.
Bonini
,
M.
Calandra
,
R.
Car
,
C.
Cavazzoni
,
D.
Ceresoli
,
G. L.
Chiarotti
,
M.
Cococcioni
,
I.
Dabo
,
A.
Dal Corso
,
S.
de Gironcoli
,
S.
Fabris
,
G.
Fratesi
,
R.
Gebauer
,
U.
Gerstmann
,
C.
Gougoussis
,
A.
Kokalj
,
M.
Lazzeri
,
L.
Martin-Samos
,
N.
Marzari
,
F.
Mauri
,
R.
Mazzarello
,
S.
Paolini
,
A.
Pasquarello
,
L.
Paulatto
,
C.
Sbraccia
,
S.
Scandolo
,
G.
Sclauzero
,
A. P.
Seitsonen
,
A.
Smogunov
,
P.
Umari
, and
R. M.
Wentzcovitch
, “
Quantum espresso: A modular and open-source software project for quantum simulations of materials
,”
J. Phys.: Condens. Matter
21
,
395502
(
2009
).
93.
Y.
Yao
and
Y.
Kanai
, “
Plane-wave pseudopotential implementation and performance of scan meta-GGA exchange-correlation functional for extended systems
,”
J. Chem. Phys.
146
,
224105
(
2017
).
94.
Y.
Shao
,
L. F.
Molnar
,
Y.
Jung
,
J.
Kussmann
,
C.
Ochsenfeld
,
S. T.
Brown
,
A. T. B.
Gilbert
,
L. V.
Slipchenko
,
S. V.
Levchenko
,
D. P.
O’Neill
,
R. A.
DiStasio
, Jr.
,
R. C.
Lochan
,
T.
Wang
,
G. J. O.
Beran
,
N. A.
Besley
,
J. M.
Herbert
,
C.
Yeh Lin
,
T.
Van Voorhis
,
S.
Hung Chien
,
A.
Sodt
,
R. P.
Steele
,
V. A.
Rassolov
,
P. E.
Maslen
,
P. P.
Korambath
,
R. D.
Adamson
,
B.
Austin
,
J.
Baker
,
E. F. C.
Byrd
,
H.
Dachsel
,
R. J.
Doerksen
,
A.
Dreuw
,
B. D.
Dunietz
,
A. D.
Dutoi
,
T. R.
Furlani
,
S. R.
Gwaltney
,
A.
Heyden
,
S.
Hirata
,
C.-P.
Hsu
,
G.
Kedziora
,
R. Z.
Khalliulin
,
P.
Klunzinger
,
A. M.
Lee
,
M. S.
Lee
,
W.
Liang
,
I.
Lotan
,
N.
Nair
,
B.
Peters
,
E. I.
Proynov
,
P. A.
Pieniazek
,
Y.
Min Rhee
,
J.
Ritchie
,
E.
Rosta
,
C.
David Sherrill
,
A. C.
Simmonett
,
J. E.
Subotnik
,
H.
Lee Woodcock
 III
,
W.
Zhang
,
A. T.
Bell
,
A. K.
Chakraborty
,
D. M.
Chipman
,
F. J.
Keil
,
A.
Warshel
,
W. J.
Hehre
,
H. F.
Schaefer
 III
,
J.
Kong
,
A. I.
Krylov
,
P. M. W.
Gill
, and
M.
Head-Gordon
, “
Advances in quantum chemical methods and algorithms in the Q-chem 3.0 program package
,”
Phys. Chem. Chem. Phys.
8
,
3172
(
2006
).
95.
N. D.
Mitri
,
S.
Monti
,
G.
Prampolini
, and
V.
Barone
, “
Absorption and emission spectra of a flexible dye in solution: A computational time-dependent approach
,”
J. Chem. Theory Comput.
9
,
4507
4516
(
2013
).
96.
C.-I.
Song
and
Y. M.
Rhee
, “
Development of force field parameters for oxyluciferin on its electronic ground and excited states
,”
Int. J. Quantum Chem.
111
,
4091
4105
(
2011
).
97.
E.
Heid
,
P. A.
Hunt
, and
C.
Schröder
, “
Evaluating excited state atomic polarizabilities of chromophores
,”
Phys. Chem. Chem. Phys.
20
,
8554
8563
(
2018
).
98.
T.
Graen
,
M.
Hoefling
, and
H.
Grubmüller
, “
AMBER-DYES: Characterization of Charge Fluctuations and Force Field Parameterization of Fluorescent Dyes for Molecular Dynamics Simulations
,”
J. Chem. Theory Comput.
10
(
12
),
5505
5512
(
2014
).
99.
G.
Tiberio
,
L.
Muccioli
,
R.
Berardi
, and
C.
Zannoni
, “
How does the trans–cis photoisomerization of azobenzene take place in organic solvents?
,”
ChemPhysChem
11
,
1018
1028
(
2010
).
100.
S.
Pipolo
,
E.
Benassi
,
G.
Brancolini
,
M.
Valášek
,
M.
Mayor
, and
S.
Corni
, “
First-principle-based md description of azobenzene molecular rods
,”
Theor. Chem. Acc.
131
,
1274
(
2012
).
101.
J. C.
Slater
and
J. G.
Kirkwood
, “
The van der waals forces in gases
,”
Phys. Rev.
37
,
682
697
(
1931
).
102.
W.
Tang
,
E.
Sanville
, and
G.
Henkelman
, “
A grid-based bader analysis algorithm without lattice bias
,”
J. Phys.: Condens. Matter
21
,
084204
(
2009
).
103.
Y. K.
Law
and
A. A.
Hassanali
, “
Role of quantum vibrations on the structural, electronic, and optical properties of 9-methylguanine
,”
J. Phys. Chem. A
119
,
10816
10827
(
2015
).
104.
M. R.
Hilaire
,
D.
Mukherjee
,
T.
Troxler
, and
F.
Gai
, “
Solvent dependence of cyanoindole fluorescence lifetime
,”
Chem. Phys. Lett.
685
,
133
138
(
2017
).
105.
H.
Lami
, “
Presence of a low-lying ‘rydberg’ band in the vapour phase absorption spectra of indole and 1-methyl indole
,”
Chem. Phys. Lett.
48
,
447
450
(
1977
).
106.
E. H.
Strickland
,
J.
Horwitz
, and
C.
Billups
, “
Near-ultraviolet absorption bands of tryptopan. Studies using indole and 3-methylindole as models
,”
Biochem.
9
,
4914
4921
(
1970
).

Supplementary Material