Absorption and Emission Spectral Shapes of a Prototype Dye in Water

Article pubs.acs.org/JPCA

Absorption and Emission Spectral Shapes of a Prototype Dye in Water by Combining Classical/Dynamical and Quantum/Static Approaches Alessio Petrone,†,#,¶ Javier Cerezo,‡,# Francisco J. Avila Ferrer,‡,§ Greta Donati,† Roberto Improta,*,∥ Nadia Rega,*,†,⊥ and Fabrizio Santoro*,‡ †

Dipartimento di Scienze Chimiche, Università di Napoli ‘Federico II’, Complesso Universitario di M.S. Angelo, via Cintia, I-80126 Napoli, Italy ‡ CNR−Consiglio Nazionale delle Ricerche, Istituto di Chimica dei Composti Organometallici (ICCOM-CNR), UOS di Pisa, Area della Ricerca, via G. Moruzzi 1, I-56124 Pisa, Italy § Physical Chemistry, Faculty of Science, University of Málaga, Málaga, 29071, Spain ∥ CNR−Consiglio Nazionale delle Ricerche, Istituto di Biostrutture e Bioimmagini (IBB-CNR), via Mezzocannone 16, I-80136 Napoli, Italy ⊥ Italian Institute of Technology, IIT@CRIB Center for Advanced Biomaterials for Healthcare, Largo Barsanti e Matteucci, I-80125 Napoli, Italy S Supporting Information *

ABSTRACT: We study the absorption and emission electronic spectra in an aqueous solution of N-methyl-6oxyquinolinium betaine (MQ), an interesting dye characterized by a large change of polarity and H-bond ability between the ground (S0) and the excited (S1) states. To that end we compare alternative approaches based either on explicit solvent models and density functional theory (DFT)/ molecular-mechanics (MM) calculations or on DFT calculations on clusters models embedded in a polarizable continuum (PCM). In the first approach (ClMD), the spectrum is computed according to the classical Franck−Condon principle, from the dispersion of the time-dependent (TD)DFT vertical transitions at selected snapshots of molecular dynamics (MD) on the initial state. In the cluster model (Qst) the spectrum is simulated by computing the quantum vibronic structure, estimating the inhomogeneous broadening from statespecific TD-DFT/PCM solvent reorganization energies. While both approaches provide absorption and emission spectral shapes in nice agreement with experiment, the Stokes shift is perfectly reproduced by Qst calculations if S0 and S1 clusters are selected on the grounds of the MD trajectory. Furthermore, Qst spectra better fit the experimental line shape, mostly in absorption. Comparison of the predictions of the two approaches is very instructive: the positions of Qst and ClMD spectra are shifted due to the different solvent models and the ClMD spectra are narrower than the Qst ones, because MD underestimates the width of the vibrational density of states of the high-frequency modes coupled to the electronic transition. On the other hand, both Qst and ClMD approaches highlight that the solvent has multiple and potentially opposite effects on the spectral width, so that the broadening due to solute−solvent vibrations and electrostatic interaction with bulk solvent is (partially) counterbalanced by a narrowing of the contribution due to the solute vibrational modes. Qst analysis evidences a pure quantum broadening effect of the spectra in water due to vibronic progressions along the solute/solvent H-bonds.

1. INTRODUCTION

challenges remain open and mainly concern the modeling of either floppy molecules or large systems with remarkable nonadiabatic couplings (in this field modern quantum-dynamical

UV−vis electronic spectroscopy is a powerful technique to investigate molecular electronic states and their interactions with the environment. Thanks to continuous progresses, computational methods have reached a sufficient maturity to help the interpretation of experimental spectra and provide a detailed understanding of the molecular characteristics that determine their position, width, and shape. However, a number of © XXXX American Chemical Society

Special Issue: Jacopo Tomasi Festschrift Received: October 31, 2014 Revised: February 20, 2015

A

DOI: 10.1021/jp510838m J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

solute−solvent explicit interactions is to resort to hybrid models, considering a cluster that comprises the solute and the directly bound solvent molecules, embedded in a polarizable continuum that describes the bulk solvent effects. However, even this approach is not straightforward. On one hand, the definition itself of the most representative cluster is, in many cases, not easy because the first-solvation sphere is not a static concept and it changes with time. On the other hand, even when a cluster model is established, it is far from being granted that a Qst approach based on a harmonic expansion of the PES can be applied and provide robust and reliable results. In fact, in principle, the PES describing the interaction of the solute and the solvent may be characterized by strong anharmonicities and multiple local minima. On these premises, here we report a detailed study of the absorption and emission electronic spectra in water of a prototype dye, N-methyl-6-oxyquinolinium betaine (MQ), characterized by a large change of polarity and H-bond ability between the ground and the excited state. We calculate the absorption and emission spectra according to the Qst and ClMD protocols and compare their predictions on the spectral shapes and the Stokes shift. To the best of our knowledge, the combination of these two approaches has been rarely pursued in literature, at least for large dyes in aqueous solution.30−32 The aim of this paper is 2-fold: first, we want to verify if and to what extent Qst calculations based on model (harmonic) PESs can be realistically applied to solute−solvent clusters chosen from MD analysis; second, we want to explore in detail which kind of information can be gained by comparing Qst and ClMD in protic solvents and if such a comparison can improve the understanding of the different roles played by the solute, the solvent, and their interactions in determining the spectral shape. To dissect the differences inherent to the different classical and quantum description of the spectra from those arising from different solvent models, we start from an analysis of the spectra in the gas phase. Later on, we proceed to the comparison of spectra in solution trying to investigate the relative role played by intramolecular (vibrational) and intermolecular (H-bond vibrations, inhomogeneous broadening) contributions to the spectral shape. Concerning the Qst calculations on cluster models based on harmonic PESs, we preliminary validate the robustness of the spectra predictions and we report a detailed analysis on the influence of the choice of the cluster model on the spectral position and shape.

approaches like MCTDH are allowing remarkable improvements1), and the accurate simulation of spectra in polar hydrogen-bonding (H-bond) solvents. In this paper we focus on the effect of H-bond solvents. To that end we consider a rather rigid molecule that can be described within the Born−Oppenheimer approximation. For this kind of system, even when they are made by several dozens of atoms, the computation of the electronic spectral shape in the gas phase is now at hand thanks to recent time-independent (TI)2−8 and time-dependent (TD)9−17 quantum approaches. These methods are quantum (Q) in the sense that they compute the manifold of the vibronic transitions associated with the change of the electronic state. This is done either by a TI sum-over-state approach or by a TD approach, Fourier-transforming the quantum correlation function between the initial thermally populated states and the corresponding doorway states propagating on the final-state potential energy surface (PES). The initial- and final-state PESs of the optical transition are usually modeled by a Taylor expansion around their equilibrium positions. In this sense, these Q methods rely on a static (st) characterization of the potential energy surfaces in terms of stationary points, gradients, and Hessians; in the following we name them Qst. The expansion in most of the cases is limited is to quadratic terms even if anharmonic effects along specific degrees of freedom can be accounted for (see, for example, refs 14 and 18−21). The adoption of implicit solvent models, like the polarizable continuum model (PCM),22 provides a straightforward route to extend the Qst approach to spectra in polar aprotic solvents. In fact, they give the possibility to account for solvent effects on the position of the spectra and their shape, considering both the modification of the underlying vibronic structure23−26 and the inhomogeneous broadening.16,27,28 When the solvent is capable of establishing hydrogen bonds with the solute, new challenges come into play and they are inherent to the necessity to account for explicit solute−solvent interactions leading, in principle, to a dramatic increase of the complexity of the solute−solvent system. The most natural approach in these cases is based on explicit solvation models, coupled with quantum mechanics/molecular mechanics (QM/MM) PESs. However, several (hundreds of) solvent molecules must be included to have a physically reliable description and, due to the complexity of this model, no straightforward application of the aforementioned Qst approach is possible. In fact, on one hand, the PES becomes very complicated with several minima and this makes unfeasible a static characterization and, on the other hand, it is not possible either to compute the vibronic eigenstates necessary for TI approach or to propagate the quantum states as required in the TD approach. Classical (Cl) molecular dynamics (MD) is the ideal approach to deal with the complexity of explicit solvent models because it is able to ensure an adequate sampling of the phase space of the PES of the initial state of the electronic transition. On these grounds it is possible to simulate the shape of the electronic spectrum by applying a classical analogue of the Franck−Condon (FC) principle,29 i.e. considering the spreading of the vertical transitions of a representative sample of snapshots of the MD trajectory. The inherent limitation of this approach, named hereafter ClMD, is that it cannot describe vibrationally resolved spectra. Moreover, although it introduces the effect of vibrational motion on the spectra, it underestimates the intrinsic delocalization of the nuclear wave function, except in the limit of high temperatures. A possible way out to save the quantum features of the electronic spectrum and account, at least, for the most important

2. THEORETICAL FRAMEWORK AND COMPUTATIONAL DETAILS In a TD formalism, the absorption and emission electronic spectra can be written as the Fourier transform of a time correlation function, χ(t,T)10,12 L(ω) =

1 2π Z v

∫ Trvi[⟨vi|⟨i|e(it−β)H /ℏμe−itH /ℏμ|i⟩|vi⟩]e±iωt dt i

f

(1)

where plus and minus signs hold for absorption and emission respectively, |i⟩ and |f⟩ are the initial and final electronic states, Tr refers to the trace operation over the vibrational states |vi⟩ associated with |i⟩, Zv is the corresponding vibrational partition function, β = (KBT)−1, KB is the Boltzmann constant, T the absolute temperature, and μ is the dipole operator written in the spectral representation μ = [μif|i⟩⟨f| + μfi|f⟩⟨i|]. The vibronic Hamiltonians associated with |i⟩ and |f⟩ are written as a sum of B

DOI: 10.1021/jp510838m J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A the relative PES (Vi or Vf) and the nuclear kinetic operator TN, Hi = TN + Vi and Hf = TN + Vf. Equation 1 is the starting expression for Qst calculations and makes evident that the quantum spectrum depends on the dynamics of the initial vibrational states on the final-state PES. Considering a set of nuclear coordinates Q, moving to a coordinate representation and neglecting TN, it is possible to get a classical approximation of the spectrum12,29 L(ω) =

1 2π

are not harmonic in principle, alternative harmonic models can be constructed depending on the nuclear structures selected to perform the quadratic expansion.34,35 More specifically, whereas the initial-state PES is expanded around its minimum, the finalstate PES is always expanded around its own minimum in adiabatic (A) approaches and around the initial-state minimum in vertical (V) approaches. The usage of the terms “adiabatic” and “vertical” in the context of vibronic spectra was first proposed in ref 36 where the underlying equations were reported. To account for Duschinsky mixings and frequency changes, it is necessary to compute the Hessian of both the initial and final PESs; for this reason in refs 37 and 38 we proposed to name these models adiabatic Hessian (AH) and vertical Hessian (VH) to help to distinguish them from simpler models where differences in normal modes and/or in their frequencies are neglected. AH and VH models have also been named respectively adiabatic Franck−Condon and vertical Franck−Condon in recent literature.18 Calculations were performed with a development version of the code -*classes .39,40 The spectra were computed in a TD approach that delivers fully converged spectra at any temperature (including 0 K) and exploits the fact that the expression for the correlation is analytical in harmonic approximation.9−16,41,42 A TI calculation of the spectrum based on a effective prescreening method4−6 was used to obtain the stick spectra of the most relevant transitions and assign them in terms of number of quanta of the different modes. Such a TI method is also implemented in Gaussian 09.8,43 Optimization and frequency calculations were performed using the Gaussian09 program package44 by means of (TD)-Density Functional Theory (DFT), adopting the PBE0 functional45 along with the 6-31+G(d,p) basis set. Solvent effects were included through the C-PCM method,22 using Pauling radii to define the cavity and, adopting, for excited states, the linear-response (LR) approach. We also considered clusters of the solute MQ and a few explicit water molecules (two, three, and four molecules around the carboxylic group of MQ) surrounded by a polarizable continuum also described by C-PCM. The solvent inhomogeneous polar broadening, due to the rearrangement of the solvent as a response to the change in the charge distribution in the solute, is accounted for through a Gaussian convolution to the vibrational structure. According to Marcus formulation,27 the standard deviation of the Gaussian (σpol) is related to the solvent reorganization energy (λ) as follows,

∫ dQ ρi (Q,T )|μif (Q)|2 ∫ dt ei(V (Q)−V (Q)±ℏω)t /ℏ i

f

(2)

where ρi(Q,T) is the density of the initial vibrational states |vi⟩ at temperature T. Noticing that the integral over the time represents a delta-function in the energy domain, we finally get L(ω) = ℏ

∫ dQ ρi (Q,T )|μif (Q)|2 δ(Vi(Q) − Vf (Q) ± ℏω) (3)

Equation 3 is the starting expression for ClMD calculation of the spectrum and it is worthy to highlight that it does not contain anymore information on the dynamics on the final-state PES. Running an MD on the initial-state PES, Vi(Q), it is possible to get a classical approximation of ρi(Q,T), ρi,Cl(Q,T). Moreover, considering a number of snapshots (Q*) that are representative of ρi,Cl(Q,T), the spectrum can be obtained from the dispersion of the transition energies |Vf(Q*) − Vi(Q*)| because, due to the presence of the δ function in eq 3, each structure Q* will only contribute to the intensity of the spectrum at the energy ℏω = |Vf(Q*) − Vi(Q*)|. The above derivation helps to highlight the differences between the Qst and ClMD approaches to the calculation of the spectra. In the quantum calculations Qst, the density ρi(Q,T) is obtained from the quantum vibrational states populated according to the Boltzmann distribution and, as noticed above, the spectrum reflects a quantum propagation on the final electronic state f. The vibronic structure of the Qst spectrum arises from these characteristics. As a drawback, for large systems as those considered in this paper, the calculations require the adoption of harmonic models for Vi and Vf PESs. In the classical ClMD approach, the vibronic resolution is lost and the density of states is approximated classically. However, because the Vi and Vf can be computed on the fly, the harmonic approximation is not needed and this allows us to introduce anharmonic effects both on the density and on the transition energies. Finally, it is worthy to notice that in harmonic models, the quantum distribution ρi(Q,T) corresponds to the spatial part of the Wigner distribution and is analytical T

ρW (Q,T ) ∝ exp[−Q CQ]

σpol 2 = 2λKBT

As reported previously,16,28 λ for absorption (λabs) and emission (λemi) can be computed as the difference of nonequilibrium (neq) and equilibrium (eq) energies by means of state-specific (SS) implementation of PCM/TD-DFT46,47 (see the Supporting Information for further details). Vertical transition energies and solvent reorganization energies were computed both using the C-PCM variant in combination with Pauling radii and using the default implementation (IEF-PCM) and default cavity in Gaussian 09. As shown in the Supporting Information (section S2), the latter protocol provides more reasonable estimates of λ and of the corresponding broadening, and therefore it was used to compute the vertical energies and inhomogeneous broadening adopted to produce the spectra compared with experiment. Actually, application of eq 6 in combination with standard Gaussian 09 SS-PCM/ TD-DFT calculations has already been shown to nicely

(4)

where C is a diagonal matrix whose elements are Ωi tanh[Ωi/ (2KBT)] and Ωi is the harmonic frequency of mode i. A simple classical approximation of the density in harmonic approximation is instead obtained through the following Boltzmann distribution ρB (Q,T ) ∝ exp[−QTΩ2Q/(2K BT)]

(6)

(5)

where Ω is the diagonal matrix of the frequencies Ωi. 2.1. Qst Approach. The quantum vibronic spectra were simulated in harmonic approximation, adopting Cartesian coordinates to describe the normal modes and accounting for differences in the normal modes and frequencies of the initial and final states (Duschinsky mixings33). Because the electronic PESs C

DOI: 10.1021/jp510838m J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

analysis reported in ref 54. Also for emission spectra, we disentangled structural and energetic solvent effects by simulating a gas phase emission spectrum by using MQ configurations extracted from the aqueous solution AIMD trajectory.

reproduce the broadening effects on the absorption spectra of a number of systems in aprotic polar solution, including coumarin C153,28 MQ,16 the chromophore of Green fluorescent protein31 and a number of additional coumarins (like c343, nkx-2311, nkx-2398, nkx-2586, nkx-2753, unpublished work). 2.2. Absorption and Emission Spectra with the ClMD Approach. To calculate absorption and emission spectra according to eq 3, ab initio molecular dynamics (AIMD) were performed for MQ in the gas phase and aqueous solution, in both the ground S0 and the excited S1 electronic states. The MQ molecule was treated at the DFT level, by using the global hybrid B3LYP functional and the 6-31G(d,p) basis set. To simulate MQ in aqueous solution, we adopted a hybrid explicit/implicit representation of the solvent, by enforcing nonperiodic boundary conditions through the GLOB model.48,49 The barycenter of the MQ solute was fixed at the center of a sphere with radius of 11 Å, including 130 water molecules (about three solvent shells) described with the standard TIP3P50 water model. Both the ground-state AIMD simulations in the vacuum and in solution were performed by exploiting the atom centered density matrix propagation (ADMP) formalism for MQ.51−53 The valence orbital’s fictitious mass μ was 0.1 amu bohr2, whereas core orbitals were weighted according to the tensorial fictitious mass scheme described in ref 52. An average temperature of 298 K was enforced by a velocity rescaling procedure. Further details about AIMD trajectories in aqueous solution are given in ref 54. The S0 AIMD was collected for 10 ps, after 4 ps of equilibration, and with 0.2 fs as time step. According to CASSCF/CASPT2 (Complete Active Space Self Consistent Field/Complete Active Space Second Order Perturbation Theory) analysis the MQ S1 excited state is well described by a HOMO → LUMO single excitation and does not show significant multireference character. Relying on that, the AIMD simulation in the excited state was performed by using energy and energy derivatives calculated at the TD-B3LYP/ 6-31G(d,p) level of theory, adopting a Born−Oppenheimer formalism.55−58 This S1 AIMD was simulated for 8 ps, after 4 ps of equilibration, and with a time step of 0.2 fs. To calculate ClMD absorption spectra, 2267 and 1666 configurations were extracted at regular time intervals from aqueous solution and gas phase S0 AIMDs, respectively. According to the AIMD (see below) the MQ oxygen is directly bound to four or fewer water molecules. On these grounds, aqueous solution configurations (MQ(H2O)130) were partitioned and treated by ONIOM with an electronic embedding, including at high level the MQ molecule and the four water molecules closest to the MQ oxygen, chosen in each frame on the basis of a oxygen−oxygen distance criterium. In analogy with the Qst treatment, vertical transition energies and dipole strengths were obtained at the TD-PBE0/6-31G+(d,p)/TIP3P/C-PCM and TD-PBE0/6-31+G(d,p) level for the aqueous solution and the gas phase, respectively. We also investigated indirect (structural) and direct (energetic) solvent effects by simulating a gas phase absorption spectrum by using MQ configurations extracted from the aqueous solution AIMD trajectory. All the simulated emission spectra were obtained in a similar manner, by extracting 1570 and 1354 configurations at regular time intervals from aqueous solution and gas phase S1 AIMDs, respectively. Concerning the emission spectrum in aqueous solution, a MQ(H2O)3 cluster was treated at high level in the ONIOM partition of the whole system, on the basis of the microsolvation

3. RESULTS This section is organized as it follows. As a first step, we compare the spectra computed in the gas phase by ClMD and Qst approaches and their differences will be discussed by analyzing the vibrational density of states (VDOS) underlying the composition of the optical spectra line shape in the two cases. For what concerns the spectra in water solution, after having verified the solidity of Qst predictions obtained on the MQ(H2O)n clusters, we shall compare the computed optical spectra with their experimental counterparts, discussing some general trends. 3.1. Spectra in the Gas Phase. Figure 1 compares the ClMD and Qst spectra computed at 298 K. Qst spectra show clear marks

Figure 1. Absorption and emission spectra in the gas phase at 298 K simulated with the Qst vibronic model (AH) and classically according to ClMD at the PBE0/6-31+G(d,p) level. MD run at the B3LYP/ 6-31G(d,p) level of theory. The vertical energies computed for emission (blue, 11950 cm−1) and absorption (red, 16990 cm−1) are also shown as thick vertical lines. Notice that no phenomenological broadening is added in Qst calculations.

of vibronic structure that, according to theory, cannot be captured by ClMD calculations. On the contrary, Table 1 shows Table 1. First Moment 41, Standard Deviation (σvib), and Maximum ωmax of the Qst and ClMD Spectra Computed in the Gas Phase at 298 K Qst ClMD

abs emi abs emi

41 (cm−1)

σvib (cm−1)

ωmax (cm−1)

16490 10500 16600 10990

1830 2210 1060 1210

15630 12330 16630 11280

that the centers of gravity 41 of the ClMD and Qst spectra are rather similar both in absorption and in emission. It is interesting to notice that ClMD calculations are able to reproduce the fact that, according to Qst, 41 is always red-shifted with respect to the vertical transition (EVabs = 16 990 cm−1 and EVemi = 11 950 cm−1) D

DOI: 10.1021/jp510838m J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 2. Distribution of dimensionless displacements of the projected MD (S0; GP) snapshots on B3LYP/6-31G(d,p) normal modes described in internal coordinates. Three representative modes, which are significantly coupled with the electronic transition, are shown: a ring bending (a), a ring breathing (b) and a CC stretching (c).

and the shift is larger for emission spectra, characterized by a longer tail extending far from the band origin. In a Qst formalism it is easy to show that these shifts do not arise from displacements of the equilibrium positions and are connected to quadratic or higher-order differences of the ground- and excited-state PESs.59 Despite the above similarities, both the absorption and emission spectra computed at the Qst level are remarkably broader than those obtained with the ClMD protocol: the standard deviations of the Qst spectra are 1830 cm−1 (abs) and 2210 cm−1 (emi) whereas they are respectively 1060 and 1210 cm−1 according to ClMD simulation. Figure 1 also shows that the Qst spectra are remarkably more asymmetric than the ClMD ones and this is a pure quantum effect due to the underlying vibronic transitions. As a consequence, the differences between the maxima and the centers of gravity of the spectra are much larger for Qst than for ClMD. However, interestingly, they have the same sign, i.e., according to both approaches the spectral maximum falls between the band origin (the 0−0 transition in Qst) and the center of gravity. 3.2. Quantum and Classical Vibrational Density of States. To better understand the features underlying the differences between ClMD and Qst spectra, in this section we focus on absorption spectra and we compare the S0 VDOS computed by the two approaches. To that end, we expressed all the nuclear configurations sampled by ClMD in terms of displacements from the B3LYP/6-31G(d,p) equilibrium geometry along the dimensionless normal modes (section S3 of the Supporting Information for details), obtaining a distribution along each mode that is compared with both the classical Boltzmann (eq 5) and the quantum Wigner (eq 4) distributions in harmonic approximation. Representative results are shown in Figure 2 for three modes respectively with frequencies 277 cm−1 (a), 590 cm−1 (b), and 1625 cm−1 (c), selected among those giving rise to the largest vibronic progressions (their dimensionless displacement during the electronic excitation is respectively 0.53, 0.56, and 1.00). In addition, Figure 3 reports the standard deviations of the Wigner, Boltzmann, and ClMD distributions for each mode, as a function of the mode frequency (the corresponding distributions are reported in section S3 of the Supporting Information, Figure S2). Data refer to S0 normal modes expressed in internal coordinates, which are better suited than Cartesian ones to deal with large displacements, as those arising from the N−CH3 torsion.60 In the frequency region of the CO and ring stretchings (1500− 1700 cm−1) the ClMD distributions are in line with the Boltzmann model (see the case of the CC stretching in Figure 2c) while they are significantly narrower than the quantum Wigner one. This difference explains why ClMD spectra are remarkably narrower

Figure 3. Standard deviations of the distribution of the displacement arising from the projection on the GS normal modes of all the MD snapshots. The standard deviation of Boltzmann and Wigner distributions at 298 K are also shown.

than Qst ones, because progressions along these modes are the main responsible for the width of the Qst spectra. On the contrary, in the low-frequency region (