Computer Simulations of the Solvatochromism of ... - ACS Publications

7704

J. Phys. Chem. B 1999, 103, 7704-7719

Computer Simulations of the Solvatochromism of Betaine-30 S. R. Mente and M. Maroncelli* The PennsylVania State UniVersity, UniVersity Park, PennsylVania 16802 ReceiVed: May 11, 1999; In Final Form: July 16, 1999

Monte Carlo simulations of the pyridinium N-phenolate dye “betaine-30” in 12 solvents (20 solvent representations) were performed in order to explore the molecular basis of the ET(30) scale of solvent polarity. Ab initio (HF/6-31G*) and semiempirical (AM1 and INDO/S) electronic structure calculations were used to determine the geometry and charge distribution of betaine-30 in its S0 and S1 states. The solvent effect on the betaine absorption spectrum was assumed to derive from electrostatic interactions between the effective charge distributions of solvent molecules and the charge shift brought about by the S0 f S1 transition. Two models for this charge shift, one obtained from INDO/S calculations and the other an idealized two-site model, were used for the spectral calculations. Good agreement between simulated and observed ∆ET shifts (ET(30) values measured relative to the nonpolar standard tetramethylsilane) was found for both charge-shift models. In water and other hydroxylic solvents, the O atom of the betaine solute was observed to form moderately strong hydrogen bonds to between one and two solvent molecules. The contribution of these specifically coordinated molecules to the ∆ET shift was found to be large, (30-60%) and comparable to experimental estimates. Additional simulations of acetonitrile and methanol in equilibrium with the S1 state of betaine-30 were used to determine reorganization energies in these solvents and to decide the extent to which the solvent response to the S0 T S1 transition conforms to linear response predictions. In both solvents, the spectral distributions observed in the S0 state simulations were ∼15% narrower than those in the S1 simulations, indicating only a relatively small departure from linear behavior. Reorganization energies were also estimated for a number of other solvents and compared to values reported in previous experimental and theoretical studies.

I. Introduction

SCHEME 1

Solvatochromic scales of solvent polarity play an important role in organic chemistry as a result of their ability to correlate the solvent dependence of equilibrium constants and rates of chemical reactions.1,2 Such scales utilize solvent-induced changes in the spectra of molecular probes to measure the relative strengths of solute-solvent interactions, i.e., relative solvent “polarities”. Polarities measured in this way are generally more useful than dielectric estimates because they are sensitive to short-range phenomena not captured in dielectric measurements, for example, the polarity of nondipolar solvents such as dioxane and benzene.3,4 As described in Reichardt’s comprehensive review,5 there are now scores of solvatochromic probes that have been used to define a large number of empirical scales of solvent polarity. In part, this proliferation of polarity scales results from the fact that different scales often disagree in their relative ranking of solvents. Such disagreement is not surprising. Although we often speak of “solvent polarity” as if it were a single, well-defined quantity, in reality there is no single solvent property that can gauge the complex array of interactions responsible for solvation. Any one solvatochromic probe merely reports “polarities” which reflect its own particular sensitivities to different solvent attributes such as electronic polarizability, dipole density, and/or hydrogen-bonding ability. For this reason, intelligent use of empirical polarity scales rests on understanding what mix of interactions is reflected in a given scale. In the present paper, we explore the nature of the interactions underlying what is probably the most popular empirical scale * Corresponding author. E-mail: [email protected].

of solvent polarities, the ET(30) scale.6,5,7 The ET(30) scale is based on the absorption spectral shifts of the pyridinium N-phenolate betaine dye “betaine-30” and its tert-butylated derivative “PTBB”, shown in Scheme 1. (PTBB is used as a surrogate for the primary probe betaine-30 in weakly polar solvents where the latter is insoluble.) Between the extremes of tetramethylsilane and water, the S0 f S1 absorption maximum of betaine-30 shifts from 930 to 453 nm, an energy change of 43.4 kcal/mol or 11000 cm-1. This enormous shift renders betaine-30 one of the most sensitive probes of solvent polarity

10.1021/jp991549r CCC: $18.00 © 1999 American Chemical Society Published on Web 08/20/1999

Simulations of the Solvatochromism of Betaine-30 known. To date, it has been used to characterize well over 300 pure solvents, a broad array of solvent mixtures and ionic solutions,5 as well as microheterogeneous media such as micelles,5 and solid and liquid surfaces.8 Most studies of betaine30 have entailed peak shift measurements for purposes of polarity characterization. However, there have also been a number of more recent studies aimed at achieving a fundamental understanding of betaine-30’s unique solvatochromism. Several experimental investigations have focused on the shape of the betaine absorption spectrum in order to derive estimates for the reorganization energy of the S0 f S1 transition.9-11 Most recently, resonance Raman spectroscopy has been combined with absorption band shape analysis in order to more accurately separate the solvent and intramolecular contributions to the spectral line shape.12,13 This interest in reorganization energies derives from recent attempts to model the S0 T S1 transition as an electron transfer process.9-19 Especially noteworthy in this regard are the studies of Barbara and co-workers,10,11,14-16 which concern the dynamics of the ultrafast internal conversion (S1 f S0) of betaine-30. Although we do not treat dynamics in the present work, we will simulate reorganization energies, an important ingredient in this dynamical modeling, and compare our results to available experimental estimates. A number of theoretical and computational studies of the solvatochromism of betaine-30 have also been reported.20-26,4,27-30 Two investigations, by Beckett and Dawber21 and de Alencastro et al.24 examined betaine-30 interacting with one or a few solvent molecules. Surprisingly good correlations were observed between ET(30) values and ground state interaction energies21 and transition frequencies.24,31 These observations provide some indication of the likely importance of specific interactions in determining ET(30) polarities. Other studies used semiempirical quantum mechanical approaches coupled to continuum23 or dipole lattice28,29 descriptions to examine the behavior in bulk solvents. The latter studies, by Lipinski and Bartkowiak, were the first to point out the important role of internal torsional coordinates in modulating the characteristics of the S0 f S1 transition. Of most interest relative to the present work are three theoretical treatments of betaine solvatochromism, which modeled betaine-30 in bulk solvents using molecular solvent models. The simplest of these treatments, reported by Richert and coworkers,26,25 employed nonpolarizable, dipolar hard-sphere representations for the solute and solvent. These workers showed that application of the MSA theory32 to such systems was able to correlate the transition energies of betaine-30 in some 110 solvents more accurately than continuum dielectric models.33 The reasonable success of this treatment would seem to imply that interactions between the permanent charge distributions of betaine and the solvent are paramount in determining the spectral shifts observed. However, a second theoretical study employing similar idealized solute/solvent models led to a rather different interpretation. Matyushov et al.27 considered how induction and dispersion interactions (not included in Richert’s treatment) affect the spectral shifts and concluded that dispersion and dipolar interactions are of nearly equal importance in determining the net shifts. The important role played by dispersion interactions in this case rests on assuming a marked increase (nearly 2-fold) in the polarizability of betaine-30 upon electronic excitation. We believe that such a large polarizability change is unlikely and inconsistent with at least some of the experimental evidence.34 Nevertheless, Matyushov et al. were able to successfully reproduce the observed transition energies of betaine-30 in a wide range of polar aprotic solvents. They noted

J. Phys. Chem. B, Vol. 103, No. 36, 1999 7705 that ET(30) values of protic solvents were not reproduced by their theory, presumably due to the importance of specific hydrogen bonding in these cases. The final statistical mechanical study reported by Perng et al.4 is quite close in spirit to the work performed here. Rather than using point dipole approximations, these authors modeled both the solute and solvent as collections of interaction sites of the sort typically employed in computer simulations. Solvation energetics were derived in several approximations using structural input from XRISM calculations.35 The solvents in their study were represented in full or nearly full atomic-level detail, but for computational expedience Perng et al. employed a simplified four-site model for betaine-30. Although one might doubt the adequacy of such an approximation, intelligent choice of interaction sites and a slight tuning of the solute charge distribution provided good agreement between calculated shifts and observed ET(30) values for nine varied polyatomic solvents. The solvents examined included both nondipolar solvents such as benzene, as well as the protic solvents methanol and water. The near quantitative agreement achieved for the variety of polar aprotic solvents studied strongly suggests that that the primary interactions responsible for the spectral shifts of betaine involve the S0 f S1 charge shift and the average charge distributions of solvent molecules. Their success in modeling methanol and water with the same treatment further suggests that nonelectrostatic aspects of hydrogen bonding are not especially important in treating interactions of betaine-30 with protic solvents. The work reported herein builds on the successes achieved by Perng et al.4 We employ Monte Carlo simulations to examine the solvation and solvatochromic shifts of a more realistic model of betaine-30 in 12 different solvents (20 solvent models). Surprisingly, with the exception of the dynamical study in acetonitrile by Lobaugh and Rossky, currently in press,30 these simulations are the first to be reported on betaine. The absence of previous simulations of this important probe undoubtedly reflects the potential difficulties involved in realistically modeling the spectroscopy of a solute as large and complex as betaine30. Among the many issues of concern are the likely inaccuracies in descriptions of the excited states of a molecule of this size, as well as the possible effects of internal flexibility and the coupling between the electronic polarizability of the molecule and the solvent environment. Several of these issues are nicely addressed in the work of Lobaugh and Rossky.30 However, the results of the Perng study suggest that the essential features of betaine-30’s solvatochromism can be understood without accounting for such effects. In the present simulations, we therefore employ a simplified representation consisting of a rigid, classical, nonpolarizable, all-atom description of betaine30 and similar models for solvent molecules. We find that such representations are indeed capable of semiquantitatively reproducing ET(30) polarities in a wide range of solvents, even in the absence of any adjustable parameters. The agreement with experimental results justifies use of these relatively unsophisticated models to interpret the molecular-level interactions responsible for betaine-30’s remarkable solvatochromy. Of particular interest is the role that specific hydrogen bonding plays in determining the spectral shifts. This interest is motivated by previous experimental studies36,37 in which we noted a good correlation between the rates of solvent-mediated proton transfer and ET(30) polarity in a large collection of aliphatic alcohols. We conjectured that, especially in alcohols, variations in ET(30) reflect mainly differences in hydrogen bond donating ability

7706 J. Phys. Chem. B, Vol. 103, No. 36, 1999

Mente and Maroncelli

rather than nonspecific “polarity”. We have therefore included eight different hydroxylic solvents in the present study. In addition, we also compare results obtained with the two different parametrizations of the alcohol solvents that were previously used in modeling these proton-transfer reactions38,39 in order to help assess the relative accuracy of these solvent models. The organization of the remainder of the paper is as follows. In section II, we discuss correlations between experimental values of ET(30) and other empirical measures of solvent polarity in order to motivate the spectroscopic model used for simulation. The intermolecular potentials and details of the simulation methodology are then briefly discussed in section III. A number of ab initio and semiempirical electronic structure calculations were performed in order to develop models for the S0 state of the solute and for the charge redistribution in the S0 f S1 transition. The results of these calculations are described in section IV. Although we employ a rigid model of betaine-30 for simulation, in this section we also examine the effects of internal torsional flexibility on the S0 and S1 states using semiempirical methods. The simulation results are then finally discussed in sections V and VI. We begin by describing the structural features and energetics of solvation of the ground state of betaine-30 and then compare the simulated spectral shifts to experimental results. In section VI, we also discuss the linearity of the S0 T S1 process and compare the reorganization energies determine from simulated line widths to previous estimates made on the basis of the experimental and theoretical models mentioned above. II. Empirical Correlations: Implications for Spectral Modeling Solvatochromic shifts originate from the same solute-solvent interactions responsible for ground state solvation. One may therefore approximately decompose these shifts in the manner40

∆ET ) hν - hν0 ) ∆Usr + ∆Uel + ∆Uind + ∆Udisp (1) where hν0 denotes the gas-phase energy difference, ∆U the difference in solvent-solute interaction energies in the two electronic states of the transition, and the designations “sr”, “el”, “ind”, and “disp” denote short-range, electrostatic, chargeinduced, and dispersion interactions, respectively. In the simulations presented here, we make the drastic assumption that the single term ∆Uel dominates the solvent dependence of the absorption shifts of betaine-30. Specifically, we assume that the difference in the transition energies of betaine in a given polar solvent and a nonpolar reference solvent, tetramethylsilane (TMS), is determined solely by the interactions between the permanent charge distributions of the solute and solvent, described by ∆Uel:

∆ET(obs) ≡ hν[polar] - hν[TMS] = ∆Uel

(2)

In so doing, we completely neglect contributions from dispersion interactions, which might arise from differences in the solute’s polarizability in the S1 and S0 states. We also ignore possible contributions that could result from changes in the effective size or shape of the solute upon electronic excitation. Induction interactions are not included explicitly either. However, by performing simulations with charge distributions parametrized for liquid-state simulations, some effects of induction interactions are included in an average sense, by virtue of the higher charges/dipole moments of such parametrizations relative to gasphase values.

Figure 1. Comparison of normalized ET(30) polarities (eq 3) to other empirical measures of solvent polarity: (top) reorganization energies of coumarin 153 (“C153 λ”); (middle) dielectric continuum reaction field factors (eq 4); (bottom) “acceptor numbers”. Filled symbols denote polar aprotic solvents and open symbols protic solvents. Open circles denote n-alcohols, squares other monohydroxy alcohols, triangles water and polyols, and hexagons other protic solvents such as amides. See Table 1 for data sources and for statistics on the regression lines plotted here.

Ultimately, the utility of the approximations embodied in eq 2 will be judged by the extent to which the simulation results reproduce experiment. Nevertheless, it is also useful to appreciate the extent to which this simple approach is motivated by the large experimental database available on betaine-30. We therefore pause to compare betaine-30 shifts to other empirical indicators of solvent polarity. For these comparisons we employ ENT values, which are normalized representations of transition energies of betaine-30, defined for a given solvent x by N

ET [x] )

hν[x] - hν[TMS] hν[H2O] - hν[TMS]

(3)

The top panel of Figure 1 compares ENT values to solvent reorganization energies (λ) measured with the polarity probe coumarin 153 (C153 hereafter). The latter quantities derive from a combined analysis of steady state and dynamic Stokes shifts data.3,42 Reorganization energies such as these are sensitive only to nuclear (as opposed to electronic) motions of solvent molecules and thus measure what we will term the “nuclear polarizability” of solvents. In principle, nuclear polarizabilities could involve dispersion interactions, but in the case of C153 λ appears to be dominated by interactions between the charge distributions of the solute and solvent molecules.3,42,43 The approximate proportionality (see Table 1) between ENT and λ(C153) in 33 aprotic solvents (filled symbols) displayed in Figure 1 therefore recommends eq 2. (The deviations exhibited by protic solvents, shown as open symbols, are discussed below.) It is worth noting that the relation of ENT with reorganization energies is not unique to C153.

Simulations of the Solvatochromism of Betaine-30

J. Phys. Chem. B, Vol. 103, No. 36, 1999 7707

TABLE 1: Regression Results Comparing ETN and Other Measures of Solvent Polaritya property

solvents excluded

N

C153 λ protic 33 reaction field protic, quadrupolar 44 acceptor number protic, halogenated 27 aprotic, acids 14

slope

R

data source

2400 cm-1 1.64 42.1 52.2

0.930 0.952 0.944 0.972

3 79 48 80

a Normalized values of ET(30) were taken from refs 5 and 7. “N” is the number of solvents, “slope” is the slope of the (zero intercept) correlation line shown in Figure 1, and R is the correlation coefficient of the fit.

Ravi et al. have recently shown that ENT values exhibit a similarly high degree of correlation (Rg 0.93 for six of the solutes) with the Stokes shifts of 10 other solvatochromic probes in selected solvents.44 Further evidence that only the nuclear polarizability of solvents is important in determining ENT is provided by the good correlation between these values and dielectric continuum predictions for the orientational component of the dipolar reaction field,45

f(0,n) )

 - 1 n2 - 1  + 2 n2 + 2

(4)

This comparison is shown in the middle panel of Figure 1. The filled symbols here denote data in some 44 aprotic solvents which contain a variety of functional groups: aldehydes, ketones, esters, ethers, nitroalkanes, and nitriles, as well as halogenated and sulfur-containing compounds. The only aprotic solvents excluded from this collection are those having small dipole moments but large higher order moments such as benzene, since dielectric continuum models are inapplicable in these cases.3 The correlation coefficient between ENT and f(0,n) for this collection of solvents is 0.952. Including a term involving solvent electronic polarizability (the n2 portion of eq 4) to the regression, does not significantly improve the fit (R ) 0.958). This insensitivity to dielectric measures of electronic polarizability, as well as the close correlation with reorganization energies noted above, leads to the conclusion that the solvatochromism of betaine-30, or equivalently the ET(30) scale of solvent “polarity”, is more or less ignorant of solvent electronic polarizability. It is useful to note that this feature of the betaine30 shifts is not universal, or even the norm. For example, similar analysis of the absorption shifts of the solute 4-nitroanisole, used to define the popular π* polarity scale,46 reveals a greater sensitivity to solvent electronic polarizability than to nuclear polarizability. (Linear regressions of the same solvent set used in Figure 1 yield correlation coefficients of 0.78 and 0.98 without and with inclusion of the electronic polarizability term.) Given the approximate proportionality between ENT and f(0,n), it is useful to examine the magnitude of the proportionality constant in light of dielectric continuum models of solvatochromism. Assuming the solute to be a point dipole centered within of a spherical cavity of radius a, continuum models predict that the slope dνjabs/df(,n) should be given by 2µ bg(µ bg b µe)/hca3, where b µg and b µe are the ground and excited state dipole moments of the solute and h and c are Planck’s constant and the speed of light. Using the experimental values b µg ) 14.8 µe) ) -8.7 D of betaine-26,47 assumed collinear, and (µ be - b and the observed slope of 6910 cm-1, one calculates an effective solute radius a ) 5.7 Å. This value is in reasonable agreement with other solvatochromic estimates of an effective spherical size for betaine-30 (for example 6.3-6.5 Å,27 5.8 Å,26 and 6.2

Å44). However, this value is significantly smaller than the estimate of 7.8 Å that is obtained from the van der Waals volume of betaine-30 (503 Å3; au ) 4.9 Å) corrected for the radius of a solvent molecule using the approximation aeff = au + jrν. (rjν denotes an average solvent radius, 2.9 Å for the set of solvents considered.) The large discrepancy between the observed and anticipated values of aeff undoubtedly reflects the difference between the extended charge distribution of betaine30 and the point dipole model assumed in deriving the dielectric continuum predictions. As will be discussed in the following section, the actual charge distributions in S0 and S1 produce much larger spectral shifts than a (solute centered) point dipole representation would predict (i.e., larger by a factor of (7.8/ 5.7)3 ) 2.6). The data in Figure 1 also illustrate the well-known fact that the shifts of betaine-30 are sensitive to the hydrogen bond donating ability of protic solvents. The ENT values of such solvents (open symbols) are significantly larger than expected based on the correlations with either f(0,n) or λ(C153) in aprotic solvents. Deviations of this sort are expected in dielectric correlations when specific solute-solvent hydrogen bonding is important, since bulk dielectric parameters are necessarily ignorant of specific solute-solvent interactions. The reorganization energies of C153 are also believed to be unaffected by specific solute-solvent interactions.42 The deviations illustrated in these two comparisons indicate that hydrogen bonding greatly enhances the shifts observed in protic solvents compared to what they would be in the absence of specific interactions. For example, using the f(0,n) correlation in aprotic solvents to calibrate the effect of nonspecific polar interactions on ENT , the portion of the shift attributable to hydrogen bonding appears to be 40-60% in most alcohol and amide solvents. Presumably, the oxygen atom of betaine-30, which is both highly charged and solvent accessible, is the primary site for this hydrogenbonding interaction. This assignment is supported by the excellent correlation observed between ENT and the “acceptor number” (AN) scale of Mayer and co-workers,48,49 shown in the bottom panel of Figure 1. The acceptor number scale is based on the 31P chemical shifts (or equivalently the PdO stretching frequencies50-52) of triethylphosphine oxide. It was originally constructed to provide a measure of solvent “electrophilicity”. Although our own analysis53 indicates that the observed NMR/IR shifts in aprotic solvents are also sensitive to nonspecific polarity effects, AN values primarily respond to hydrogen bonding interactions between the solvent and with the O atom of the probe. The fact that a good correlation with ENT is found in aprotic solvents and an excellent correlation is found in protic solvents (see Table 1) indicates that both of these polarity scales register a blend of nonspecific and specific solute-solvent interactions, and that the mix is only slightly different in the two cases.45 This similarity underscores the importance of specific interactions at the O site of betaine-30 in determining its UV shifts. III. Simulation Methods and Potentials The simulations reported here are Monte Carlo calculations carried out using the “BOSS” molecular simulation program developed by Jorgensen.54 Simulations were performed in the isothermal-isobaric (NPT) ensemble at a temperature of 25 °C and pressure of 1 atm. The systems studied consisted of a single solute molecule immersed in a cubic box of solvent molecules. The box size (32-34 Å on a side) was chosen such that a 10 Å layer of solvent would surround the solute on each side after the system had reached equilibrium. The number of molecules

7708 J. Phys. Chem. B, Vol. 103, No. 36, 1999

Mente and Maroncelli

required ranged from 357 for tert-butyl alcohol to 1249 molecules for water. Cubic periodic boundary conditions were applied and solvent-solvent interactions were spherically truncated at a cutoff distance based on oxygen-oxygen atom distances. A solute-solvent cutoff was applied such that if any solute atom/solvent oxygen atom distance was less than this cutoff value the interactions between the entire solute and solvent molecule were included. The values of both of these cutoffs were chosen to be between 12 and 13 Å in all of the simulations. No attempt was made to account for electrical interactions beyond these cutoff distances. It is likely that proper incorporation of long-range interactions would slightly alter the magnitudes of the solvation energies simulated for the highly dipolar betaine solute. However, by maintaining the same, relatively large system size for all solvents, variations among different solvents, which are our primary interest, should be little affected. All simulations were started from an initial random configuration that had been equilibrated from a liquidlike arrangement of solvent molecules for a period of at least 108 configurations. Energies and densities were monitored in order to ensure adequate convergence within the equilibrium period. Following equilibration, simulations were performed in sets of five or more segments of 2 × 107 configurations each, to compute statistical uncertainties. The uncertainty values reported here are (1 times the standard deviation of the mean of the series of segments. Molecules were represented as collections of interaction sites and intermolecular interactions modeled via site-site LennardJones plus Coulomb terms:

[( ) ( ) ]

νij ) 4ij

σij rij

12

-

σij rij

6

+

qiqj rij

TABLE 2: Solute and Selecteda Solvent Potential Parameters species

no.

betaine-30

source 55 56

CCl4

1

54

chloroform

2

81

acetone

3

54

acetonitrile

5o

82

acetonitrile

5r

83

water (TIP4P) 12o

84

atom or group

q i/ au

σii/ Å

ii/kB/ kcal mol-1

C N H(-C) O(-C) C Cl CH Cl O C CH3 C N CH3 C N CH3 O H M

b b b b 0.248 -0.062 0.420 -0.14 -0.424 0.300 0.062 0.28 -0.43 0.15 0.129 -0.398 0.269 0.000 0.520 -1.040

3.55 3.25 2.42 2.96 3.80 3.47 3.80 3.47 2.96 3.75 3.91 3.650 3.200 3.775 3.40 3.20 3.60 3.1536 0.0000 0.0000

0.07 0.17 0.03 0.21 0.050 0.266 0.08 0.30 0.210 0.105 0.160 0.150 0.170 0.207 0.0993 0.170 0.380 0.155 0.000 0.000

a Parameters for all solvents not shown in the above table have been previously reported in ref 39. b See part B of Table 3.

(5)

The Lennard-Jones parameters between unlike atoms were derived from the like-atom parameters (provide below) using the mixing rules σij ) (σiiσjj)1/2 and ij ) (iijj)1/2. In general, each atom within a molecule corresponds to an individual site, with the main exception being the CHn groups of the alcohol solvents, which were taken as single units centered on the carbon atom. The betaine-30 solute was represented using a rigid all-atom model. Its Lennard-Jones parameters were taken from the OPLS potential functions for aromatic rings55 or nucleotide bases.56 They are summarized in Table 2 for convenience. The geometry employed was that optimized for the ground state of betaine30 with the semiempirical AM1 Hamiltonian, using the AMPAC program.57 Charges for the ground state were determined from electrostatic potential fits58 to the ab initio HF/6-31G* wave function, using the Gaussian-94 program suite.59 (More discussion of these charges is provided in the following section.) Two different solvent representations were used in this work. Both representations are based on the OPLS force field54 and employ the bond lengths, bond angles, torsional potentials, and Lennard-Jones parameters of the standard OPLS set. They differ only in the choice of site charges. In one representation, referred to hereafter as the “OPLS” solvent models, the generic charges of the OPLS liquid-state parametrizations were employed. In the second representation, charges were derived from electrostatic potential fits to ab initio (HF 6-31 G*) calculations of the isolated solvent molecules. Solvents modeled using these calculated charges are “ab initio” solvents. The motivations behind the use of these two charge representations as well as values of the Lennard-Jones and charge parameters for most of the solvents studied were previously provided in refs 38 and

Figure 2. Space-filling representation of the betaine-30 (AM1) geometry used in simulation. The darkest shading shows atoms N6 and O17 (see Figure 3 for numbering).

39 and will not be reproduced here. Parameters of those solvents not described in these references are summarized in Table 2. IV. Electronic Structure of Betaine-30 The geometry of the betaine molecule employed in the simulations is the minimum energy structure of the AM1 Hamiltonian illustrated in Figure 2. Of note is the high degree of nonplanarity of the molecule, which results from twisting of the five phenyl rings pendant to the pyridinium N-phenolate core as well as twisting about the central torsional angle of this core “χ” (see Figure 3). As will be discussed shortly, the χ angle is of particular relevance for the spectroscopy of betaine-30. In the (AM1) minimum energy conformation χ ) 48° (relative to zero at coplanar pyridinium and phenolate rings). This value

Simulations of the Solvatochromism of Betaine-30

J. Phys. Chem. B, Vol. 103, No. 36, 1999 7709

Figure 3. Schematic depiction of the ground state charges and S1-S0 charge difference used in simulation. The ground state charges (“q0”) are atomic charges fit to reproduce the electrostatic potential of the ab initio wave function. The charges labeled “∆q” are the charge difference predicted by semiempirical calculations (see text). The charges are proportional to the volumes of the spheres shown on each atom. Black represents negative and white positive charge. (The structures of betaine shown here have been flattened in order to better display the charges.)

may be compared to the experimental value of 65° observed in the crystal structure of a Br-substituted derivative of betaine30.60 Other recent semiempirical calculations by Bartowiak and Lipinski29 and Lobaugh and Rossky30 provided values of χ ≈ 60° and χ ) 56°, respectively. Some characteristics of the ground state charge distribution of betaine-30, derived from ab initio calculations (HF/6-31G*// AM1), are listed in Table 3A. The calculated dipole moment is 16.4 D, 10% larger than the experimental value of the related molecule betaine-26.47 (Semiempirical calculations by the two groups mentioned above yielded values of 16.828 and 15.8 D.30) Table 3A also lists the Mulliken charges determined from the

ab initio calculations. These charges do not exhibit the O-/N+ character anticipated from simple valence bond arguments. Instead, both the O and N atoms have large excess Mulliken populations. However, the surrounding C atoms in the “p1” and “p2” rings (see Figure 3) are also highly charged, and render a net negative and positive character to the phenolate (“p1 + O” in Table 3A) and pyridinium (p2) rings. The semiempirical calculations are similar in this respect. For simulation purposes, we employ atom-centered charges chosen to best represent the electrostatic potential of the ground state ab initio wave function using the RESP method.58 These charges are depicted schematically in Figure 3 (“q0”), and a

7710 J. Phys. Chem. B, Vol. 103, No. 36, 1999

Mente and Maroncelli

TABLE 3: Calculated Charge Distributions of the S0 and S1 States (A) Various Electrical Propertiesa state

source 6-31G* INDO/S AM1 exptb INDO/S AM1/CI exptb INDO/S AM1/CI exptb

S0

S1 S1-S0

ET/kcal mol-1

33.1 35.2 27.0 (864 nm) (809 nm) (1055 nm)

MT/D

µ/D

qO/au

qN/au

qp1/au

qp2/au

qp1+O/au

-0.71 -0.72 -0.37

-0.77 -0.13 +0.12

+0.40 +0.29 -0.14

+0.05 +0.22 +0.30

-0.31 -0.43 -0.51

-0.61 -0.40

-0.33 -0.04

-0.44 -0.09

+0.79 +0.41

-0.51 -0.60

7.6 3.5 3-6

16.4 18.5 12.9 14.8 3.9 14.9 6.2 -14.6 +3.7 -8.7

0.11 -0.04

-0.20 -0.15

-0.73 -0.12

+0.57 +0.17

-0.62 -0.16

(B) Charge Differences Employed in the Spectral Calculationsc (S1-S0) charge differences no.

atom

S0 6-31G*

INDO/S

two-site

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

C C C C C N H H C C C C C C H H O

0.0099 -0.1396 0.0032 -0.1396 0.0100 0.2349 0.1590 0.1590 -0.0802 -0.1918 -0.2324 0.5736 -0.2322 -0.1919 0.1772 0.1773 -0.6100

-0.1538 -0.0018 -0.1527 -0.0016 -0.1555 -0.1831 0.0000 0.0000 +0.3099 -0.0506 +0.1163 +0.0349 +0.1167 +0.0069 0.0000 0.0000 +0.1144

0 0 0 0 0 -0.327 0 0 0 0 0 0 0 0 0 0 +0.327

a All calculated properties refer to the AM1 minimum energy geometry. E is the transition energy, M the transition dipole, and µ the permanent T T dipole moment. The charges (“q”) listed here are all Mulliken charges. “p1” and “p2” denote the phenolate and pyridinium rings as defined in Figure 3. b Experimental values are from ref 5. The energy is the estimated gas-phase transition energy of betaine-30. Dipole moment values are from electrochromic measurements of the related molecule betaine-26 in dioxane solution.47 The experimental transition moments decrease with increasing ET and have values in the range of 3-6 D.64 c Atom numbering as defined in Figure 3. Charges listed for S0 are RESP-fit charges, whereas the charge differences are from Mulliken charges. All charges are in atomic units.

partial listing is provided in Table 3B. The RESP-fit charges differ from the Mulliken charges in that positive charge is assigned to the N atom. The most notable feature of the latter charge distribution is the large charge separation ((0.6e) in the C12-O17 bond, which produces a local dipole moment of 3.5 D. It is this dipole moment, coupled to the solvent through the accessible O atom, that gives rise to strong specific interactions between betaine-30 and protic solvents. As discussed in section II, we assume that the solvent dependence of the betaine absorption spectrum reflects only the excitation-induced changes to the electrostatic component of the solvent-solute interactions. We model these changes using differences between the effective atomic charges of the S1 and S0 states. Unfortunately, determination of excited state charges (and other properties) is not as straightforward as in the case of the ground state. Whereas the S0 charge distribution can be reliably calculated using ab initio methods and a good basis set, comparable calculations for the S1 state require a much greater investment of computational resources and faith. As an alternative, we estimate the requisite S1-S0 charge differences using configuration interaction calculations (“CI”) with the semiempirical AM161 and INDO/S62 Hamiltonians. The AM1CI calculations were performed on 396 configurations (singly and doubly excited) selected from 14 MOs using the automatic configuration selection of the AMPAC program.57 The INDO/S calculations were performed using the ARGUS program.63 The CI in this case involved 225 singly excited configurations

selected from 30 MOs. Both calculations were performed at the AM1 ground state geometry. Characteristics of the S0 and S1 charge distributions obtained from these semiempirical calculations are also provided in Table 3 and illustrated in Figure 3. The accuracy of the transition information provided by such calculations can be partially judged by comparing transition energies (ET) and transition moments (MT) to experimental values (Table 3A). One finds that both methods overestimate the transition energy by about 4 kcal/mol (∼3000 cm-1) compared to the estimated gas-phase value. The transition moments measured in the experiment range between 3 and 6 D in most solvents, with values that generally decrease with increasing ET.64 Both calculations yield values close to this range, with the INDO/S value perhaps being in better agreement with what would be extrapolated for the gas phase. In view of the uncertainties in the experimental values, these transition properties can be said to be reasonably well predicted by both semiempirical methods. However, the same is not true of the electrical properties of S1. The AM1/CI calculations predict an increase in dipole moment between the S0 and S1 states, which would lead to small solvatochromic shifts of the opposite sign from those observed in experiment. The charge redistribution predicted by this method (Figure 3) must therefore be viewed as unsuitable for our purposes and will not be considered further. On the other hand, the INDO/S dipole moments are qualitatively in line with what is known from experiment. The S0 dipole moment is slightly larger and

Simulations of the Solvatochromism of Betaine-30

J. Phys. Chem. B, Vol. 103, No. 36, 1999 7711 TABLE 4: Parameters Characterizing the Internal Torsion Angle Dependence of Various Propertiesa property

a0

a2

a4

a6

est. error

S0 energy S1 energy S0 dipole moment S1 dipole moment S0 f S1 transition energy S0 f S1 transition moment

5.176 29.140 18.762 3.374 34.315

2.396 8.273 -2.565 2.176 10.669

4.655 -4.847 0.455 -0.860 -0.193

1.841 -2.780 0.258 -0.304 -0.939

0.900 1.456 0.175 0.263 1.245

7.142

5.704

-0.884

0.289

0.627

a Parametrized fits of the properties obtained from INDO/S calculations as a function of torsional angle (see text). Units are kcal/mol for energy quantities and Debye for dipole and transition moments. a0-a6 are the coefficients in the Fourier series of eq 6 for χ in degrees. “Est. error” is the standard error of the fit.

Figure 4. Dependence of various properties on the internal torsional angle (χ) as predicted by INDO/S calculations. The curves in this figure are fits to eq 6. Values of the parameters are provided in Table 4. The dark symbols mark the angle of the AM1 ground state minimum energy geometry used in simulation.

the S1 value slightly smaller than experimental estimates, such that the large decrease in dipole moment between S0 and S1 is captured, albeit in exaggerated form. Thus, the INDO/S ∆q distribution shown in Figure 3 should provide a reasonable starting point for modeling. The charge shift entails movement of roughly 0.6e from the phenolate (“p1 + O”) to the pyridyl (“p2”) ring (Table 3A). As shown in Figure 3, this charge shift is not localized to the O and N atoms alone but rather it is distributed over four important atoms in each ring. We therefore employ the INDO/S Mulliken charge differences shown in Figure 3 to model the solvent effect on the transition energy. For the sake of simplicity, we include only charges within the core structure (atoms 1-17); the actual charge differences used are provided in Table 3B.65 Since these INDO/S charges are not expected to be quantitatively accurate, we also explore the spectral shifts simulated using a second, simplified model for the charge shift in the S0 f S1 process. In this alternative, “two-site” model, we assume that the charge shift takes place only between the O and N heteroatoms of betaine and choose the amount of charge (0.33e) to reproduce the experimentally observed dipole moment change of 8.7 D. Comparisons between these two shift models provide an indication of the effect that inaccuracies in the INDO/S charge differences might have on the simulation results. A final aspect of these gas-phase calculations concerns the influence of the central torsional angle χ (Figure 3) on transition properties. To explore the likely magnitude of such effects, we have repeated the INDO/S calculations at a series of AM1 geometries optimized for fixed values of χ. The results of these calculations are presented in Figure 4. The curves in Figure 4 are fits of the calculated properties to an expansion of the form

A(χ) ) a0 + a2 cos(2χ) + a4 cos(4χ) + a6 cos(6χ) (6) Values of the fitting parameters are provided in Table 4. Not surprisingly, the minimum energy geometry (χ ) 40°) of the S0 state predicted by the INDO/S calculations differs slightly

from the AM1 minimum (χ ) 48°, solid symbols). Between the limits of χ ) 0° and 90° one finds that the energies of the S0 and S1 states vary nonmonotonically, but in such a manner that the transition energy decreases continuously with χ from 43.8 to 24.4 kcal/mol. The permanent dipole moment of the S0 state increases 25% between χ ) 0° (µ ) 16.9 D) and χ ) 90° (µ ) 21.5 D). The S1 state dipole varies in a contrary manner, such that the average of the two moments is nearly independent of angle. Finally, the S0 T S1 transition moment is predicted to decrease monotonically with increasing χ to a value that is nearly zero at χ ) 90°. Several of the above results can be compared to calculations recently reported by Bartkowiak and Lipinski,29 who used a slightly different INDO/S method. We note that the transition energies, dipole moment changes, and oscillator strengths obtained from the present calculations agree with the gas-phase results in Figure 3 of this earlier work to 20% or better. One important observation concerning the χ dependence illustrated in Figure 4 is that the ground-state torsional potential is sufficiently shallow that there should be a nonnegligible spread of torsional angles present at room temperature. Considering only the χ coordinate, in thermal equilibrium at 300 K the rms spread in torsional angle is calculated to be (5°. This relatively narrow distribution translates into a narrow distribution of ground state dipole moments (or dipole differences) of only (0.4 D. However, the effect on the absorption spectrum is significant. The absorption spectrum (∝P(χ) ET(χ) |M(χ)|2) is predicted to be approximately Gaussian in shape with a width of σ ) 1.9 kcal/mol (fwhm ) 1600 cm-1) or a reorganization energy (λ ) σ2/2kBT, see section V) of 1100 cm-1. Although not included in the present simulations, this intramolecular source of broadening would be expected to add to any solventinduced broadening of the absorption spectrum.30 Finally, it should also be noted that the above considerations all apply to betaine-30 in the gas phase. One might expect some variation of these “internal” parameters in solution. For example, the increase in the dipole moment of S0 with χ, coupled to the relatively small difference in gas-phase energy between the minimum of the S0 potential and χ ) 90° might drive the system to a 90° geometry in polar solvents. Such behavior was indeed predicted by the Langevin dipole modeling of betaine-30 in water reported by Bartkowiak and Lipinski.29 These authors found χ ≈ 90° to be most stable in solution. This change in χ, together with the solvent’s interaction with the solute electronic polarizability, leads to an enhancement of the predicted S0 dipole moment from the value of 17 D in the gas phase to 24 D in “water”. A remarkably similar change was also found by Lobaugh and Rossky.30 In coupled classical/semiempirical simulations of betaine-30, the latter authors found an enhanced

7712 J. Phys. Chem. B, Vol. 103, No. 36, 1999

Mente and Maroncelli

TABLE 5: Characteristics of Betaine Solvationa solvent

EtN

no.

model

CCl4 chloroform acetone tert-butyl alcohol

0.052 0.259 0.355 0.389

acetonitrile

0.460

2-propanol

0.546

1-propanol

0.617

ethanol

0.654

methanol

0.762

ethylene glycol TFEb water

0.790 0.898 1.000

1 2 3 4a 4o 5o 5r 6a 6o 7a 7o 8a 8o 9a 9o 10 11 12a 12o 12r

OPLS OPLS OPLS ab inito OPLS OPLS EMM ab inito OPLS ab inito OPLS ab inito OPLS ab inito OPLS ab inito ab inito ab inito TIP4P TIP3P

acetonitrile methanol

5o 9o

OPLS OPLS

-VUV

-VUV(LJ)

% coul

72 89 96 85 83 126 149 108 99 113 100 95 96 103 106 134 93 119 128 127

72 74 67 58 60 74 82 61 64 72 60 54 56 60 60 83 59 48 47 47

0 17 30 32 28 41 45 44 36 36 40 43 42 42 43 38 37 60 63 63

1.7 1.0

11 11

21 13

1.4 1.4 1.2 1.8 1.7 1.4 1.4 1.5 1.0 1.5 2.0 2.1 2.2

13 11 11 11 8 10 10 10 13 13 7 8 8

13 16 10 17 13 18 19 14 10 20 14 13 14

S1 State Simulationsd 107 -0.70 +0.43 80

74 59

31 27

1.0

5

6

q(O)

q(H)

-0.77 -0.70

+0.43 +0.43

-0.76 -0.70 -0.72 -0.70 -0.71 -0.70 -0.67 -0.70 -0.67 -0.61 -0.79

+0.43 +0.43 +0.42 +0.43 +0.42 +0.43 +0.42 +0.43 +0.42 +0.43 +0.39

-0.83

+0.42

NHB

-VHB

%HB

(8)c

a

Columns labeled “no.” and “solvent model” designate solvent potential functions as either “ab initio (a)” or “OPLS (o)” representations, which differ only in the atomic charges. In the case of hydroxylic solvents the OH charges (“q(O)”, “q(H)”, in atomic units) are listed to illustrate the magnitude of the differences involved. VUV denotes the solute-solvent interaction energy. VUV(LJ) is the contribution of Lennard-Jones interactions to this solvation energy and “% coul” is the percentage contributed by Coulombic interactions. NHB and VHB refer to solvent molecules hydrogen bonded to the O17 site of betaine. NHB is the average number of H-bonded solvent molecules, determined by integrating the solute-solvent O-O site-site radial distribution functions until their first minimum. VHB is the average solute-solvent pair energy of solvent molecules hydrogen bonded to the oxygen site of the solute. Finally, “% HB” denotes the fraction of the total solvation energy (VUV) contributed by these H-bonded molecules. All energies are in units of kcal/mol. Statistical uncertainties (1σ level) in these values are approximately (5% for the total energy quantities (VUV) and (0.3 for the coordination numbers and hydrogen-bonding energies. b TFE denotes 2,2,2-trifluoroethanol. c In the case of chloroform, which is represented by a united (CH) site in the OPLS representation, no clear peak in the pair energy distribution is observed. The value listed indicates the position of a prominent shoulder in the distribution. d Results from simulations in which the charges of the betaine solute were taken to be the ab initio S0 values plus the INDO/S (S1 - S0) charge differences listed in Table 3.

S0 dipole moment of ∼25 D in acetonitrile as compared to 16 D. (Lobaugh and Rossky did not observe a significant change in χ as a result of solvation.) In light of these two results, it is reasonable to expect that the S0 dipole of betaine-30 is significantly larger in polar solvents than the present gas-phase calculations indicate. This fact should be kept in mind when viewing the Monte Carlo simulations described below, which are all run using the ab initio gas-phase charge distribution. V. Simulations of Solvation and Solvatochromism Before considering the spectral shifts, we first examine some energetic and structural features of the ground state solvation of betaine-30. Table 5 summarizes the results obtained for the set of 12 solvents studied here. This solvent set consists of four polar aprotic and eight hydroxylic solvents. In the case of the latter liquids, we compare the results obtained with two slightly different potential models, labeled “OPLS” and “ab initio”. As described in section III (and ref 39), the models differ only in their atomic charges, the most important of which are provided in Table 5. The column labeled “VUV” in Table 5 lists the average solutesolvent interaction energies. These solvation energies lie in the range 70-130 kcal/mol, unusually large for a nonionic solute. Such energies are mainly the result of the number of solutesolvent contacts made possible by the large size of the betaine molecule. For most of the solvents studied here one estimates that approximately 30 solvent molecules comprise the first solvation shell of betaine-30. The number is much larger in the case of water, ∼65. In solvents of moderate to high polarity (as judged by ENT ), the total solvation energy consists of

approximately equal contributions from Lennard-Jones and Coulombic interactions (“% coul” in Table 5). Although there is no clear relationship between the total interaction energy and the value of ENT , there is a modest correlation between this polarity measure and the Coulombic contribution to the energy. Given the structure of betaine-30 discussed in section III, one anticipates that the primary site for specific solvent interactions will be the phenolic oxygen atom. This expectation is verified by radial distribution functions such as those illustrated in Figure 5. Hydrogen bonds are made to the O17 site in all of the protic solvents examined. (Strong correlations such as these are not found for other solute sites, for example N6). The peak in all O-O distribution functions occurs at 2.70 ( 0.03 Å, in good agreement with the value of 2.71 Å observed in the crystalline ethanol complex of a closely related betaine.60 Measures of the number and strength of the hydrogen bonds made to the O17 site are provided in the columns labeled “NHB” and “VHB” in Table 5. The values of NHB are coordination numbers obtained by integrating the first peak of the solute-solvent O-O radial distribution functions to the first minimum position (∼3.5 Å). The values of VHB are the average pair interaction energies between the solute and a hydrogen-bonded molecule (defined as any molecule counted in NHB). Both the O-O distances, and the values of VHB ≈ 11 kcal/mol are indicative of hydrogen bonds of moderate strength.66 In most of the protic solvents either one or two solvent molecules are coordinated to O17 at any instant such that on average the coordination number is 1.4 ( 0.3. Exceptions are water, ethylene glycol, and tert-butyl alcohol. The small size of water allows for an average coordination number of 2, independent of solvent representation.

Simulations of the Solvatochromism of Betaine-30

Figure 5. Representative radial distribution functions observed in the OPLS representations of tert-butyl alcohol, methanol, and water. The solid curves are solute-O17/solvent-H atom distribution functions and the dashed curves are solute-O/solvent-O functions.

In ethylene glycol, only a single hydrogen bond is made. Finally, in the case of tert-butyl alcohol the predicted coordination numbers differ markedly between the two solvent representations. In the ab initio representation NHB is larger than in other alcohols, whereas in the OPLS representation it is smaller (despite the similar values of VHB). This difference appears to be real and must derive from the difference in O charges in the two solvent models. Apparently, the larger oxygen charge in the ab initio solvent is sufficient to overcome the steric difficulty of coordinating two bulky tert-butyl alcohol molecules to the betaine O-terminus (see Figure 2). The final column in Table 5 (“% HB”) provides an estimate of the contribution of specifically coordinated solvent molecules to the total solvation energy, using the simple product NHB × VHB. One finds that H-bonded solvent molecules account for only 15 ( 4% of the total solvation energy of the betaine molecule. As will be shown presently, these molecules contribute to a much greater extent in determining the spectral shifts. We now consider what the simulations predict for the spectra of betaine-30. Figure 6 (top panel) illustrates simulated “absorption spectra” in several solvents using the INDO/S model for the S0 f S1 charge redistribution. These spectra are simply distributions of the solvent-induced shift in the energy gap ∆ET ) Uel(S1) - Uel(S0) observed during the course of a simulation. Such distributions represent the solvent contribution to the absorption spectrum in the inhomogeneous limit, (σ∆E/ p)τ∆E . 1, where σ∆E and τ∆E represent the width and characteristic time for equilibration of the ∆ET distribution. Many prior studies of solvation dynamics67-69 have shown that the fastest components of solvation have time constants greater than 100 fs in all solvents except water (where τ∆E > 20 fs). Even with these limiting fastest times together with the smallest values of σ∆E observed, one finds uncertainty products greater than 3. It is therefore safe to assume the validity of the inhomogeneous limit and interpret these shift distributions directly as the solvent contribution to the absorption spectrum. To characterize these distributions, both the first (〈∆ET〉) and second

J. Phys. Chem. B, Vol. 103, No. 36, 1999 7713

Figure 6. Representative spectra calculated with the INDO/S chargeshift model. The top panel shows the spectra observed in a variety of solvents, designated according to the numbering scheme in Tables 5 and 6. The bottom panel shows spectra calculated in equilibrium with the ground and excited states of betaine-30 in acetonitrile (circles) and methanol (squares), along with fits of these data to Gaussian functions. Except for the ethylene glycol data (top panel, no. 10) all of the data displayed here are from simluations employing OPLS solvent representations.

TABLE 6: Characteristics of Simulated Spectral Shift Distributionsa solvent

expt ∆ET

CCl4 1.7 chloroform 8.4 acetone 11.5 tert-butyl alcohol 13.2 acetonitrile

14.9

2-propanolb

17.9

1-propanolb

20.0

ethanol

21.2

methanol

24.7

ethylene glycol TFEc water

25.6 28.8 32.4

no.

solvent model

INDO/S 〈∆E〉 σ∆E

two-site 〈∆E〉 σ∆E

1 2 3 4a 4o 5o 5r 6a 6o 7a 7o 8a 8o 9a 9o 10 11 12a 12o 12r

OPLS OPLS ab initio OPLS OPLS EMM ab initio OPLS ab initio OPLS ab initio OPLS ab initio OPLS ab initio ab initio ab initio TIP4P TIP3P

-0.5 1.2 -0.9 1.3 7.7 1.8 6.7 1.7 14.4 3.0 12.1 2.3 14.6 2.1 16.5 2.2 9.8 1.7 10.8 1.3 17.3 3.9 15.5 3.1 19.1 3.6 15.1 2.7 19.8 (4.0) 20.4 (4.1) 15.6 (3.7) 15.6 (3.6) 16.5 3.5 16.1 3.3 15.9 (3.7) 16.8 (3.7) 19.0 3.4 18.2 3.4 18.2 3.0 18.8 3.3 20.1 3.3 19.5 3.4 21.1 3.4 20.6 3.4 21.4 2.9 18.7 2.5 13.9 2.6 16.1 3.1 24.4 3.9 23.4 3.3 27.5 4.2 26.6 3.8 27.3 4.2 25.5 3.4

a Experimental values are the shifts of the absorption maximum of betaine-30 relative to the spectrum in tetramethysilane (from ref 5). “INDO/S” and “two-site” refer to shifts calculated on the basis of the charge redistribution models described in section IV. Uncertainties (1σ level) in the averages 〈∆E〉 and widths σ∆E are typically