J. Phys. Chem. B 2007, 111, 3953-3959
3953
Molecular Simulation of Solvent-Induced Stokes Shift in Absorption/Emission Spectra of Organic Chromophores Ekaterina A. Nikitina,*,† Alexey V. Odinokov,† Fedor V. Grigoriev,‡ Mikhail V. Basilevsky,† Artem A. Khlebunov,† Vyacheslav A. Sazhnikov,† and Mikhail V. Alfimov† Photochemistry Center, Russian Academy of Sciences, Moscow, Russia, and Research Computing Center, Moscow State UniVersity, Moscow, Russia ReceiVed: NoVember 16, 2006; In Final Form: January 19, 2007
The values of steady-state solvatochromic Stokes shifts (SS) in absorption/emission electronic spectra of organic chromophores are studied theoretically in the framework of the Hush-Marcus model. Charge distributions for chromophore solutes in their S0 and S1 states are found by means of conventional quantumchemical methods combined with the continuum PCM approach for treating solvation effects. The solvent reorganization energies, which are expected to correlate with the solvent-induced part of 1/2 SS, are found in a molecular dynamics (MD) simulation which invokes a novel method for separation of the inertial piece of the electrostatic response (Vener, et al. J. Phys. Chem. B 2006, 110, 14950). Computations, performed in two solvents (acetonitrile and benzene), consider three organic dyes: coumarin 153 as a benchmark system and two other chromophores, for which experimental spectra are also reported. The results are found to be in reasonable agreement with the experiment. A consistent treatment of nonlinear effect in the solvent response, promoted by the polarizability of solutes and contributing to the solvent reorganization energies (Ingrosso, et al. J. Phys. Chem. B 2005, 109, 3553), improves the results of computations.
1. Introduction Experimental measurements of time-dependent Stokes shifts in electronic emission spectra provide rich information on the relaxation dynamics of molecular solvation shells and have been intensively developing during the past decade.1-3 The theoretical interpretation, mainly focused on the relaxation phenomenon, has been suggested both in the framework of continuum solvent models4-8 or in simplified molecular theories9-12 and also based on molecular simulations.3,13-17 The attempts of theoretical estimates of absolute magnitudes of the solvent-induced component of the total steady-state Stokes frequency shift (relative to the frequency of the corresponding absorption transition) have also been reported. They are customarily based on the HushMarcus model18 and include either fitting the parameters of a semiempirical molecular-level scheme19 or computations in terms of a continuum solvent model20 as well as those based on molecular simulations.14-16 Quantum-chemical computations of pertaining charge distributions in the ground and excited solute states involved in a given spectroscopic transition were performed with different degrees of sophistication.14-16,19-21 Altogether, the list of computational studies addressed to a static solvation component of the Stokes shift is rather short and the perspectives of such calculations are not quite clear so far for several reasons. In the present work, we report and discuss the new results of molecular dynamics (MD) simulations for four organic chromophore solutes I-IV. The first one (C153) is a well-known benchmark system extensively studied both experimentally and theoretically.14,15,19,20,22 Computations for the two other species seem to be new, as well * Corresponding author. Phone:
[email protected]. † Russian Academy of Sciences. ‡ Moscow State University.
+7
095
9350120.
E-mail:
as the experimental measurements, which are therefore also briefly reported in section 6. It should be noted that the experimentally measured counterpart of chromophore III is 9-ditolylamino-2,7-dimethylacridine (IV), in which two p-CH3 groups are inserted in phenyl rings and two extra CH3 groups are present in the acridine structure. Several computations for this system have also been performed in our work and compared with similar results for system III. Internal rotations in both structures III and IV are sterically hindered. Obvious from a general consideration, this conclusion was supported by our computations. Thereby, chromophores III and IV can indeed be considered as rigid systems. MD simulations have been established as good means for modeling the temporal evolution of the Stokes shift (SS). Within such limited objectives, the problem of separation of slow (nuclear) and fast (electronic) solvent motions is not of vital importance. The situation becomes different, however, provided one is focused on absolute values of a SS or its solvation component. This is why our present MD simulation is based on a new combined continuum/molecular approach, which allowed for explicitly withdrawing the contribution of the highfrequency (inertialess) electronic component of the dielectric response from the total effect,23 in order to reveal its lowfrequency (inertial) component, which can only be observed in time-resolved spectroscopic experiments.1-3 This methodology proved to work successfully in studies of solvent reorganization energies of several electron transfer (ET) reactions.23,24 2. Basic Relations The Stokes shift (SS) is defined as a difference
SS ) hν01 - hν10
10.1021/jp067610r CCC: $37.00 © 2007 American Chemical Society Published on Web 03/27/2007
(1)
3954 J. Phys. Chem. B, Vol. 111, No. 15, 2007
Nikitina et al. The approximation λ0s ) λ1s is often applied.18-20 Being a consequence of the linear response approach, it becomes true within conventional continuum models of a solvent. Its violation at a molecular-level consideration is often not essential. A similar approximation (λ0v ) λ1v) for vibrational modes is less reliable. We restrict our study to the solvent-induced part of the SS and compare its experimentally measured values with solvent reorganization energies (λs) found in MD simulations. Thereby, the quantity to be considered is
1 λs ) (λ0s + λ1s) 2
(3)
Note that eq 2 implies that vibrational and solvent modes are separable. This is an approximation legitimate for a special class of chromophore systems in which intramolecular geometry changes during S0 f S1 transition are not large and, most importantly, they do not promote significant changes in electron charge distributions. We called such systems “rigid”, and structures I-IV are considered as rigid below. In such a case, the vibrational component in eq 2 can be treated as solventindependent for a solute in a given electronic state. It can be estimated experimentally as SS for a nonpolar aliphatic hydrocarbon solvent, in which solvent reorganization energies are negligibly small. We shall continue the further discussion in the frame of the equation
1 λs ) (SS - SS0) 2
(4)
Here, SS0 is the Stokes shift measured in such inert hydrocarbon solvent. The right-hand part of this equation defines the “experimental” solvent reorganization energy. For an adequate comparison with eq 3, average (first moment) frequencies are used in eq 1.18,22 This assumption implies also that inertialess electronic polarization effects do not contribute to eq 4, which is indeed so for the experimental Stokes shifts, owing to the fact that the kinetics of equilibration of fast electronic degrees of freedom is unobservable experimentally. On the other hand, a separation of inertialess polarization effects in the theoretically determined λs is performed in a special procedure described below. Because it is λs, which proves to be the basic quantity actually provided by the present theoretical computations, the discussion below is performed in terms of this quantity rather than in terms of the SS itself, measured experimentally. 3. Simulation of Reorganization Energies Solvent reorganization energies were computed at a molecular level as23 Frequencies ν01 and ν10 correspond to vertical transitions S0 f S1 (absorption) and S1 f S0 (emission), respectively. Potential energy surfaces for S0 and S118-20 are represented in terms of a set of nuclear coordinates, both intramolecular vibrational modes of a chromophore solute and those associated with motions of solvent particles. The initial spectral state (i.e., S0 for absorption and S1 for emission) is thermally equilibrated with respect to these coordinates. Under this condition, SS ) λ0 + λ1, where λ0 and λ1 are total reorganization energies corresponding to states S0 and S1. By dividing them into intramolecular (vibrational, index v) and outer sphere (solvent, index s) parts, we obtain
SS ) (λ0s + λ1s) + (λ0v + λ1v)
(2)
λs ) -
1 (〈F0|Φ0〉 + 〈F1|Φ1〉 - 〈F0|Φ1〉 - 〈F1|Φ0〉) 2∞
(5)
where Fi (i ) 0,1) are charge distributions for electronic states 0 (for S0) and 1 (for S1) and Φi are the corresponding solvent electrostatic response fields. Both Fi and Φi are vectors in the representation based on point atomic sites of a solute and 〈...〉 mean scalar products of these quantities, that is, the integrals in the real space. Fields Φi(r) (r is a spatial coordinate) are averages performed over equilibrium ensembles which are represented by a MD trajectory equilibrated with respect to either F0 or F1. This is a standard MD simulation of solvent structure in the presence of a solute, which does not specially account for electronic polarization of solvent particles. The scaling factor
Molecular Simulation of Stokes Shift 1/∞ in eq 5 allows for a consistent extraction of the inertial solvent response contribution from this conventional and easily available procedure. This basic equation defines λs as the inertial component of the total reorganization energy, from which its fast electronic inertialess component is eliminated. Such separation of fast and slow components of the medium polarization response field is performed in terms of a combined molecular/ continuum model,24 in which inertialess effects are considered within the dielectric continuum approximation using optical dielectric constant ) ∞, whereas inertial polarization is considered explicitly at a molecular level. As discussed earlier,24 simple and transparent formula 5 follows from this model under the assumption that long-range Coulomb interactions are treated in a uniform way in the whole space. This requirement could not be obeyed in earlier works,24 where the reaction field approach was implemented in the MD simulation. This is why a more complicated computational scheme was then applied. After applying the Ewald technique,23 eq 5 has become quite legitimate and it is now invoked in the present study. 4. Charge Distributions Charge distributions for states S0 and S1 were calculated in two steps. First, the equilibrium geometries were found. These computations, performed without solvent, used the fast and efficient density functional technique implemented in the PRIRODA package.25 The standard DFT/B3LYP/aug-cc-pVTZ method was used for S0, and it was substituted for timedependent density functional theory (TDDFT) in the same realization for S1. At the second step, charge distributions F0 and F1 were computed. Now, solvation effects were taken into account within a continuum PCM (Polarizable Continuum Model) scheme26 implemented in the Gaussian 03 (G03) package.27 We performed RHF/6-31G**/PCM computations for S0 and CIS/6-31G**/PCM computations for S1. Fixed gas phase geometries, determined individually for S0 and S1 at the first step, were used. Charge distributions F0 and F1, being transformed into ESP (Electrostatic Potential) point charge representation, were inserted in the MD simulations of section 3. The rigid solute geometry used in the simulations was determined as the average of optimized S0 and S1 gas phase geometries. In this way, the average inherent to the experimental estimate (eq 4) was simulated. A semiempirical quantum-chemical methodology was also used at the second step of charge computations. It was based on the AM1 Hamiltonian28 combined with the FRCM continuum solvation scheme.29 The FRCM scheme represents a two-cavity model which refines the standard PCM approach by considering the intermediate solvent layer in the close vicinity of the solute cavity with ) 1 inside this cavity. This intermediate layer forms a second cavity with ) ∞, whereas ) 0 (the static dielectric permittivity, as in PCM) is accepted outside this twocavity construction. The width of the second cavity is considered as an empirical parameter, a property of the solvent (it was taken as δ ) 1.8 Å29 for acetonitrile (AN)). The quantum-chemical part of semiempirical computations was based on the RHF method for the S0 state and on the self-consistent CIS/PCM scheme for the S1 state, in which the solvent response field obtained within the FRCM method was inserted. The CI expansion used (CIS 10/20) provided a satisfactory convergence of the results.
J. Phys. Chem. B, Vol. 111, No. 15, 2007 3955 5. The Nonlinearity Correction In the earlier work (see ref 19 and references therein), the importance of taking into account the solute polarizability, especially for excited states, has been emphasized. The PCM procedure, involved in a computation of the solute charges as described above, is based on the nonlinear SCRF methodology. The Schroedinger equation for the solute includes solvent response fields (Φi) in the Hamiltonian. Therefore, the charge distributions F0 and F1 are Φ-dependent; that is, the polarizability effects are introduced. However, the higher order nonlinearity still remains lost; it arises due to differences in the solute polarization for the cases of equilibrium and nonequilibrium solvation shells. For the S1 state, equilibrium (F1eq) and nonequilibrium (F1neq) distributions can be computed, which are equilibrated to inertial components of Φ(1) and Φ(0), respectively.26 Nonequilibrium charge (F1neq) corresponds to the vertical Frank-Condon 0 f 1 transition. Its computation is described in ref 15. The similar refinement for F0 charges is unavailable now due to the limitation of the pertaining G03 procedure (a HF nonequilibrium computation for S0 is impossible). In computations using eq 5, we considered only equilibrium charges. This is a conventional assumption, involved for calculations of ET reorganization energies. It is essentially based on the linear response approximation, which implies parabolic free energy profiles and a single unique charge distribution for each state. A refinement of this approach was implemented in ref 15.30 The corresponding modification of eq 5 reads
λs )
1 (〈∆F01|Φ0〉 + 〈∆F10|Φ1〉) 2∞
(6)
∆F01 ) F1neq - F0eq ∆F10 ) F0eq - F1eq In totally linear approximation, Fieq ) Fineq and eq 5 is recovered. The nonlinear correction for ∆F10 (F0eq to be changed for F0neq) is not introduced in the second term due to a technical reason mentioned above. It is expected, however, to be of secondary importance owing to the lower polarizability of the S0 state. It should also be noted that the vertical excitation producing ∆F01 must be calculated on the basis of the nonequilibrium (relative to S0) geometry,15 thus assuming that the S1 geometry has gotten relaxed and only the solvation shell remains nonequilibrium. 6. Results a. Experimental Measurements. Structure IV (9-ditolylamino-2,7-dimethylacridine) was synthesized in the Photochemistry Center, RAS.31 Nile Red dye (structure II), as well as solvents for spectral measurements (hexane, benzene, acetonitrile), were received from Aldrich; all of these products were used without further purification. The solution concentrations were 10-4 mol/L for measurements of absorption spectra and 10-5 mol/L for fluorescence spectra. The experiments were performed with a Hitachi 3100 spectrophotometer and a Hitachi 850 spectrofluorimeter. The results are displayed in Figures 1 and 2. Earlier spectral studies for Nile Red have also been reported.32 For structure IV, the spectra, both absorption and emission, have also been measured in a wider set of solvents.31 On the basis of these data and using a conventional solvatochromic
3956 J. Phys. Chem. B, Vol. 111, No. 15, 2007
Figure 1. Spectra of 9-ditolylamino-2,7-dimethylacridine: (a) adsorption; (b) fluorescence. The solvents are hexane (1), benzene (2), and acetonitrile (3).
model with the Mataga-Lippert polarity function,4 the dipole moment and its change were estimated (µ1 ) 12 D, ∆µ ) 10 D). Similar estimates for Nile Red are32a µ1 ) 14 D and ∆µ ) 9 D. Hexane was used as inert solvent, providing reference data to be inserted into eq 4. The solvatochromic shifts in absorption spectra proved to be practically solvent-independent for structure IV (Figure 1a), whereas a notable dependence does appear for the case of Nile Red (Figure 2a). The shifts in fluorescence spectra are strongly dependent on the nature of the solvent in both cases (Figures 1b and 2b). Average (first moment) frequencies corresponding to the S1 states were used for further estimations of Stokes shifts. As compared to similar estimates based on spectral maximum positions, the present results are reduced by 800-1000 cm-1. This produces the corresponding correction in solvent reorganization energies which is 2 times as small. b. Computations. The procedures described above were applied for a calculation of charge distributions F0 and F1 (QM part) and solvent reorganization energies (λs) (MD part) for four chromophore solutes I-IV. These solutes were chosen because their experimental absorption and fluorescence spectra in different solvents are available (see refs 22, 31, and 32 and section 6a) and values of solvent reorganization energies can be calculated on the basis of the experimental data. Experimental values and values calculated in MD runs of reorganization energies for two solvents, namely, acetonitrile and benzene, for all solutes under study are presented in Table 1. Both simple
Nikitina et al.
Figure 2. Spectra of Nile Red: (a) adsorption; (b) fluorescence. The solvents are hexane (1), benzene (2), and acetonitrile (3).
linear response results (eq 5) and those including nonlinear correction (eq 6) are reported for calculated λs. The calculated state dipole moments (µ0 and µ1) and absolute values of their vector differences (∆µ) are also listed. As it is seen from this table, the calculated values of reorganization energy reasonably fit experimental ones for C153 and for Nile Red for both solvents. Nonlinear correction improves the results of computations. As concerns the third molecule (9-DPAA), no agreement was obtained between experimental and calculated reorganization energies with ab initio CIS/PCM charge distributions. The underlying reason is a week degree of charge transfer for S1 found in a computation and the small value of the corresponding change of dipole moment (∆µ) (see Table 2). If one takes the more strongly polarized second excited state S2 instead of S1 and works with its charge distribution, the agreement between experimental and calculated λ’s for AN becomes better but still not satisfactory (λs ) 0.81 × 103 cm-1). Therefore, the charge distributions for 9-DPAA were received at the semiempirical AM1/FRCM level. In this case, the calculated reorganization energies of acetonitrile and benzene become closer to the experimental values. It is relevant to note that in the case of the other two solutes (structures I and II) the application of the semiempirical approach for a calculation of charge distributions did not provide good fitting to the experiment. On the basis of semiempirical charge distributions, dipole moments µ0 and µ1 and their differences were also calculated for S0 and S1. For 9-DPAA, these values are quite reasonable, whereas the ab initio calculations provide values of µ1 which are certainly too small.
Molecular Simulation of Stokes Shift
J. Phys. Chem. B, Vol. 111, No. 15, 2007 3957
TABLE 1: Solvation Part of Reorganization Energies and Solute Charge Distributions solutea
solvent
µ0, D
coumarin 153, G03
AN BZ AN BZ AN BZ AN
9.3 8.2 9.3 8.0 2.5 1.6 2.6
Nile Red, G03 9-DPAAe structure IVe
µ1 (eq/neq),b D
∆µ (eq/neq),b D
17.1/14.4 14.4 18.9/18.3 15.3 16.1 15.4 16.8
7.9/5.3 6.4 9.8/9.1 7.4 13.6 13.8 14.2
λMD (eq 5/eq 6),c cm-1 × 103
λexp ) SS/2,d cm-1 × 103
0.70/0.88 0.21 0.74/0.82 0.18 1.50 0.66 1.93
1.16 0.35 0.75 0.41 (1.57)f (0.74)f 1.57
a
Calculations based on gas phase DFT (S0) and TDDFT (S1) optimized geometries. b Based on solute charges for the S1 state computed for equilibrium and nonequilibrium inertial solvent configuration within the PCM scheme.15 c Calculation without (eq 5) and with (eq 6) polarization correction. d One-half of the solvation contribution to the SS. e Semiempirical charge calculations by AM1/FRCM (see the text). f Experimental SS values measured for structure IV.
TABLE 2: Ab initio CIS/PCM Dipole Moments and Their Changes for 9-Diphenylamino-acridine (Structure III) and 9-Ditolylamino-2,7-dimethylacridine (Structure IV) in Acetonitrilea state
excitation energy, eV
µ0, D
µex,b D
∆µ, D
structure III, S1 structure III, S2 structure IV, S1
0.5 1.6 0.4
3.1 3.1 3.4
2.1 12.3 4.3
1.3 9.2 2.1
a
TABLE 4: Computed Solvent Reorganization Energies for ET in the Cation of Structure V in Several Solventsa solvent b
water acetonitrilec dichloroethaneb benzenec sc-CO2b
G03 computations in AN. b Dipole moments for states S1 and S2. b
TABLE 3: Ab initio Gas Phase TDDFT Dipole Moments and Their Changes for 9-DPAAa µex, D
excitation energy, eV
∆µ, D
TDDFT functional
µ0, D
S1
S2
S1
S2
S1
S2
PBE0 B3LYP PBE96 BHHLYP SVWN5
0.9 1.2 1.0 1.3 1.1
17.5 18.0 18.4 16.9 18.5
0.6 0.8 2.3 1.1 3.0
16.9 17.1 17.6 15.8 17.6
0.5 0.5 1.5 0.3 2.1
2.1 1.9 1.3 2.8 1.3
3.1 3.0 2.7 3.3 2.7
a
λs, kcal/mol 25 23 17 8 5
a Note that λ = 2-3 kcal/mol for the SS of C153 in AN (Table 1). s Reference 23. c The present work.
the experiment. The position of the change transfer S1 level is more or less correct in this scheme. The main error arises, however, because the sequence of these states on the energy axis becomes inverted. Finally, in order to illustrate how eq 5 works when the ET reorganization energies are considered (a more simple case than the SS computation), we calculated values of λs for hole transfer in a cation of structure.
See ref 33.
As is seen from Table 1, the values of ∆µ closely correlate with calculated values of λ. Similar studies were performed for structure IV. Let us recall that it is this system for which SS values, using experimental estimates for λs and addressed to 9-DPAA in Table 1, were actually measured. One observes from Table 2 that CIS/PCM ab initio charge distributions are invalid for IV as well as for 9-DPAA. Computations of the reorganization energy (λs) in AN with AM1/FRCM charges confirm the close similarity of structures III (the model) and IV (the original) and of their spectral properties. In Table 3, the gas phase ab initio TDDFT computations of charge distributions for ground and excited states of 9-DPAA are represented in terms of dipole moments and their changes. The basis used was 6-31G** for all atoms except for N atoms, where 6-31G**+ was used. The dipole moments were computed using full charge distributions, contrary to the results of Tables 1 and 2, where ESP point charges were implemented. For a number of functionals, similar and stable results are obtained.33 They confirm large values of µ and ∆µ for S1. The lowest singlet excitation proves to be a charge transfer state in all cases. The irregularity in excitation energies is more significant. The dipole moments are similar to those obtained by the AM1/FRCM method (Table 1) in benzene (BZ) solvent, thus certifying the validity of this semiempirical approach, at least qualitatively. It is the locally excited state S2 for which the transition energy, as computed by the CIS/PCM method, proved to be too low. This is seen from a comparison of Table 2 with Table 3 and
These results are included in Table 4. 7. Discussion One concludes from Table 1 that the methodology based on eqs 5 and 6 gives the values of solvent reorganization energies in a satisfactory qualitative agreement with the experiment. Relative discrepancies are not small; however, the smallness of the observed effect should be borne in mind as well. The absolute deviation, which is always less than 1 kcal/mol (∼350 cm-1), is an acceptable level of accuracy for a complicated computational scheme involving inevitable elements of semiempirical reasoning. It is expedient, in this respect, to show the results of similar calculations for an ET reaction, where the reorganization energies are significantly larger. For this purpose, we list in Table 4 values of λs for the ET in the cation of structure V, calculated by means of eq 5 for different solvents. Nonlinear correction does not seem necessary for the ET case, when low-polarizable ground-state type electronic states are involved. The experimental λs value for cation V (the hole transfer) is known34 only for the reaction in dichloroethane, and its accuracy is problematic. By introducing the energy scale correction, the estimate λs = 15 kcal/mol was recommended.23 At such
3958 J. Phys. Chem. B, Vol. 111, No. 15, 2007 background, the observed deviations of 2-3 kcal/mol between the calculated and experimental values again look reasonable. On the other hand, the trends in the solvent dependence of λ’s look also quite natural, such as a clear correlation between their values and solvent dielectric constants (for polar solvents) and a similar correlation with solvent quadrupole moments22,35 (for nonpolar solvents like BZ and CO2). Returning back to the spectral-based solvent reorganization energies, we observe from Table 1 that absolute λs values constitute 2-4 kcal/mol in AN, reducing to 1-2 kcal/mol or less for the case of BZ. The ambiguity brought about by specific nuances of MD simulation techniques, including parametrization of force fields, the approximations inherent to the procedure (eq 5) performing a separation of the inertial response field, and so forth, is also expected to be of the order of 1 kcal/mol or so. Note that a typical inaccuracy of experimental measurements for the SS in the case of C153 is 200 cm-1,22 which reduces to 100 cm-1 for λs. We expect the experimental errors for other chromophores in Table 1 to be larger, mainly due to a simplified identification of the experimental SS as a shift in the spectral-averaged frequency value, and roughly estimate their upper bound again as 1 kcal/mol or 350 cm-1. The impact of the “rigid solute” approximation is also worth discussion. As noted in the Introduction, this model implies that there is no correlation between the excitation-induced changes in geometry and in charge distributions of solute chromophores. This property is strongly structure-dependent, and the relative importance of its violation has been revealed during our studies of the nonlinear correction (eq 6). In the present work, we used the gas phase geometry optimization procedure for both ground (S0) and excited (S1) states (this is consistent with the “rigidity” assumption). Let us denote the corresponding equilibrium geometries as G0 and G1. The nonequilibrium vertical S1 charge distribution, entering eq 6, was computed on the basis of G1 geometry, as noted at the end of section 5. However, provided that the rigid approximation is acceptable, the same result would appear on the basis of the alternative G0 geometry. The actual differences in µ1 and ∆µ values, based on G0 and G1 geometries and observed for the case of C153 in AN, were found to be small, but they became significant for Nile Red in AN. We illustrate this observation here by listing the pertaining changes in λs values: λs ) 840 cm-1 for the nonequilibrium C153/AN system with G0 geometry (compared to λs ) 880 cm-1 in Table 1, where G1 geometry was used) and λs ) 1110 cm-1 for Nile Red/AN (compared to λs ) 820 cm-1 in Table 1). From this example, one can conclude that Nile Red is nonrigid in larger degree than C153 and the corresponding nonrigidity effect constitutes about 300 cm-1 (versus 40 cm-1 for C153). For the BZ solvent, the geometry changes from G1 to G0 also produced slight changes in values of λs, showing the impact of the nonrigidity effect (70 cm-1 for Nile Red). Note also that the nonlinearity effect, arising owing to eq 6, completely vanishes in BZ, because the nonequilibrium shift of the S1 charge distribution, appearing due to a change of the inertial solvent polarization, is absent for the case of a nonpolar solvent (0 ) ∞), within the continuum PCM approach implemented in our G03 computations of solute charges. We can conclude that, within the present semiempirical treatment, there exists inevitably a lot of uncontrollable factors violating its validity and a rough estimate of the total ambiguity introduced in this way may gain ∼1 kcal/mol, at least for the present systems under investigation. At this background noise, regular trends, which can be seen from Table 1, point to a
Nikitina et al. reasonable correlation between calculated and experimental values of solvent reorganization energies. Because the deviations between the experiment and a computation change systematically in Table 1, they can be further reduced by a proper parametrization, but this would hardly bring new physical information, because a lot of approximations, even not mentioned here, are inherent to any theoretical scheme of such a sort. The accuracy of 1 kcal/mol in a quantum-mechanical energy computation for excited states of large organic molecules is hardly available at the CIS level. The errors of 0.5 eV for energy levels are usual.36-38 The strategy accepted in the present work is based on the underlying hypothesis that the accuracy of charge computations within the same technique may be relatively higher. Then, by using purely electrostatic eq 5 or 6, one can gain acceptable results. Thereby, the most delicate and laborious ingredient of a computational job becomes the work with charge distributions for excited states. Serious problems still exist on this way. The breakdown of the ab initio CIS method for 9-DPAA and its counterpart IV, found in the present work, seems to be typical. This is illustrated by Table 2, where unrealistically small values of dipole moments and their changes are shown for S1 and S2 states of these solutes. It seems that our hypothesis formulated above works only until the CIS scheme provides a correct order of energy levels for lowest excited states in solution. The inversion of energy levels results from exaggerated stabilization of locally excited states, which seems to be typical for extended aromatic structures. When the energy levels are inverted, the calculated polarization of the charge transfer state seems to become reduced. On the other hand, the dipole moments for C153 in Table 1 are close to those of ref 15 and agree well with the experimental estimates.39 For Nile Red and structure IV, the computed µ1 and ∆µ values are notably higher than those estimated independently by means of a conventional empirical procedure4,32 (see section 6a), but a qualitative agreement remains and a relative success of the present λs computations argues in favor of the quantum-chemical results. Special test computations of charge distributions for ground and excited states of 9-DPAA in the gas phase, performed by means of the TDDFT technique (see Table 3), confirmed the suggestion that the first exited singlet of this chromophore is a charge transfer state; this conclusion can be extended to structure IV. Altogether, one may conclude that the difference in our treatment of structures I and II, on one hand, and III and IV, on the other hand, appears at the level of a quantum-chemical computation and is a result of failure of a simplified CIS/PCM approach as applied to the latter pair of structures. In this respect, a perspective of self-consistent TDDFT/PCM computations of solute charges, when this becomes available,40 looks promising. Acknowledgment. The present work was performed with the financial support of RFBR grants 04-03-32445 and 05-0308162-ofr. The authors are greatly indebted to Prof. Benedetta Mennucci, Prof. Branka M. Ladanyi, and Prof. Marshall D. Newton for fruitful discussions. We are grateful to Dr. A.A. Granovsky for performing gas phase TDDFT computations for 9-DPAA. Appendix: Simulation Details MD trajectories were monitored in a cubic cell with constant volume (NVT ensemble) at T ) 300 K. The box length was 4 nm, and periodic boundary conditions were applied. The rigid
Molecular Simulation of Stokes Shift solute molecule was placed at the cell center, its location and orientation being fixed during the MD run. The number of solvent particles provided the experimental density value, that is, 726 molecules of acetonitrile and 431 benzene molecules, taking into account the solute excluded volume. The coupling constant with the Berendsen thermostat was 0.3 ps. The trajectory step size was 1 fs, and the typical length of the trajectory for performing averages was 200 ps with a preliminary equilibration for 100 ps. Several test computations of averages were made with extended trajectories with the length 400 ps. The corresponding changes in reorganization energies were always less than 5 × 10-3 kcal/mol. The GROMACS package41 was used in all computations. The Ewald PME technique for treating electrostatic interactions was used; the cutoff radius for Lennard Jones (LJ) potentials was 1.3 nm. Force field parametrization for solutes was borrowed from the GROMOS-96 manual42 corrected recently by Oostenbrink et al.43 LJ parameters and effective atomic site charges for benzene were extracted from the same source. The rigid nonpolarizable threecenter model of acetonitrile44 was used. References and Notes (1) Maroncelli, M. J. Mol. Liq. 1993, 57, 1. (2) Stratt, R. M.; Maroncelli, M. J. Phys. Chem. 1996, 100, 12981. (3) Ladanyi, B. M. In Theoretical Methods in Condensed Phase Chemistry; Schwartz, S. D., Ed.; Kluwer: Dordrecht, The Netherlands, 2000. (4) Bakshiev, N. G. Opt. Spektrosk. 1964, 16, 446. Mazurenko, Y. T. Opt. Spektrosk. 1974, 36, 283; 1980, 48, 388. Kakitani, T.; Mataga, N. J. Phys. Chem. 1987, 91, 6277. (5) Bagchi, B.; Oxtoby, D. W.; Fleming, G. R. Chem. Phys. 1984, 86, 257. (6) van der Zwan, G.; Hynes, J. T. J. Phys. Chem. 1985, 89, 4181. (7) Song, X.; Chandler, D.; Marcus, R. A. J. Phys. Chem. 1996, 100, 11954. Hsu, C. P.; Song, X.; Marcus, R. A. J. Phys. Chem. B 1997, 101, 2546. Song, X.; Chandler, D. J. Chem. Phys. 1998, 108, 2954. (8) Parsons, D. F.; Vener, M. V.; Basilevsky, M. V. J. Phys. Chem. A 1999, 103, 1171. (9) Wolynes, P. G. J. Chem. Phys. 1987, 86, 5133. (10) Rips, I.; Klafter, J.; Jortner, J. J. Chem. Phys. 1988, 88, 3246; 1988, 89, 4288; J. Phys. Chem. 1990, 94, 8557. (11) Bagchi, B.; Chandra, A. J. Chem. Phys. 1989, 90, 7338; 1989, 91, 2594; 1992, 97, 5126. Chandra, A.; Bagchi, B. J. Chem. Phys. 1993, 99, 553. (12) Fried, L. E.; Mukamel, M. J. Chem. Phys. 1990, 93, 932. (13) Carter, E. A.; Hynes, J. T. J. Chem. Phys. 1991, 94, 5961. (14) Kumar, P. V.; Maronchelli, M. J. Chem. Phys. 1995, 103, 3038. (15) Ingrosso, F.; Ladanyi, B. M.; Mennucci, B.; Elola, M. D.; Tomasi, J. J. Phys. Chem. B 2005, 109, 3553. (16) Ingrosso, F.; Ladanyi, B. M.; Mennucci, B.; Scalmani, G. J. Phys. Chem. B 2006, 110, 4953. Ingrosso, F.; Ladanyi, B. M. J. Phys. Chem. B 2006, 110, 10120. (17) A full list of the literature devoted to MD simulations of the timedependent Stokes shift can be found in refs 3 and 13-15. (18) Hush, N. Prog. Inorg. Chem. 1967, 8, 391. Dogonadze, R. R.; Itskovich, E. M.; Kuznetsov, A. M.; Vorotyntsev, M. A. J. Phys. Chem. 1975, 79, 2827. Marcus, R. A. J. Phys. Chem. 1989, 93, 3078. (19) Matyushov, D. V.; Newton, M. D. J. Phys. Chem. A 2001, 105, 8516. (20) Mertz, E. L. J. Phys. Chem. A 2005, 109, 44. (21) For the literature on earlier quantum-chemical computations of the solute charge distributions, see refs 15, 19, and 20. (22) Horng, M. L.; Gardecki, J. A.; Papazyan, A.; Maroncelli, M. J. Phys. Chem. 1995, 99, 17311. Reynolds, L.; Gardecki, J. A.; Frankland, S.
J. Phys. Chem. B, Vol. 111, No. 15, 2007 3959 J. V.; Horng, M. L.; Maroncelli, M. J. Phys. Chem. 1996, 100, 10337. (23) Vener, M. V.; Tovmash, A. V.; Rostov, I. V.; Basilevsky, M. V. J. Phys. Chem. B 2006, 110, 14950. (24) Leontyev, I. V.; Vener, M. V.; Rostov, I. V.; Basilevsky, M. V.; Newton, M. D. J. Chem. Phys. 2003, 119, 8024. Leontyev, I. V.; Tovmash, A. V.; Vener, M. V.; Rostov, I. V.; Basilevsky, M. V. Chem. Phys. 2005, 319, 4. (25) Laikov, D. N. Chem. Phys. Lett. 1997, 282, 151. (26) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. ReV. 2005, 105, 2999. (27) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision B.03; Gaussian Inc.: Pittsburgh, PA, 2003. (28) Stewart, J. J. P. Mopac ’97; Fujitsu Limited: Tokyo, Japan, 1999. (29) Basilevsky, M. V.; Rostov, I. V.; Newton, M. D. Chem. Phys. 1998, 232, 189. Newton, M. D.; Rostov, I. V.; Basilevsky, M. V. Chem. Phys. 1998, 232, 201. (30) The authors are grateful to Prof. B. M. Ladanyi for providing this important information. (31) Sazhnikov, V. A.; Khlebunov, A. A.; Alfimov, M. V. High Energy Chem. (in Russian) 2007, 41, 1. (32) (a) Golini, C. M.; Williams, B. W.; Foresman, J. B. J. Fluoresc. 1998, 8, 395. (b) Boldrini, B.; Cavali, E.; Painelli, A.; Terenziani, F. J. Phys. Chem. A 2002, 106, 6286. (c) Moog, R. S.; Kim, D. D.; Oberle, J. J.; Ostrovski, S. G. J. Phys. Chem. A 2004, 108, 9294. (33) The TDDFT computations for 9-DPAA have been made by Dr. A. A. Granovsky with the PCGAMESS package. (34) Johnson, M. D.; Miller, J. R.; Green, N. S.; Closs, G. L. J. Phys. Chem. 1989, 93, 1173. (35) Zimmt, M. B.; Waldeck, D. H. J. Phys. Chem. A 2003, 107, 3580. (36) Dreuw, A.; Head-Gordon, M. Chem. ReV. 2005, 105, 4009. (37) Cammi, R.; Corni, S.; Mennucci, B.; Tomasi, J. J. Chem. Phys 2005, 122, 104513. (38) Lappe, J.; Cave, R. J.; Newton, M. D.; Rostov, I. V. J. Phys. Chem. B 2005, 109, 6610. (39) Moylan, C. R. J. Phys. Chem. 1994, 98, 13513. Smirnov, S. N.; Braun, C. L. ReV. Sci. Instrum. 1998, 69, 2875. Chowdhury, A.; Locknar, S. A.; Premvardhan, L. L.; Peteanu, L. A. J. Phys. Chem. A 1999, 103, 9614. Samanta, A.; Fessenden, R. W. J. Phys. Chem. A 2000, 104, 8577. (40) Improta, R.; Barone, V.; Newton, M. D. ChemPhysChem 2006, 7, 1211. (41) van der Spoel, D.; van Bruuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A. L. T. M.; Hess, B.; Feenstra, K. A.; Linhdahl, E.; van Drunen, R.; Berendsen, H. J. C. Gromacs User Manual, version 3.1.1; Groningen, The Netherlands, 1991-2002. (42) van Gunsteren, W. F.; Billeter, S. R.; Eising, A. A. Hunenberger, P. H.; Kruger, P.; Mark, A. E.; Scott, W. R. P.; Tironi, J. G. Biomolecular Simulation: The GROMOS96 Manual and User Guide; vdf Hochschulverlag AG and der ETH Zurich and BIOMOS b.v.: Zurich, Switzerland, and Groningen, The Netherlands, 1996. (43) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. J. J. Comput. Chem. 2004, 25, 1656. (44) Edwards, D. M. F.; Madden, P. A.; McDonald, I. R. Mol. Phys. 1984, 51, 1141.
