J. Phys. Chern. 1995, 99, 955-964
955
Reaction Dynamics of Charge-Transfer State Formation of 4-(N,N-Dimethylamino)benzonitrilein a Methanol Solution: Theoretical Analyses Shigehiko Hayashi, Koji AndoJ and Shigeki Kato* Department of Chernistly, Faculty of Science, Kyoto University, Kitashirakawa, Sakyo-ku, Kyoto 606,Japan Received: May 20, 1994; In Final Form: October 25, 1994@
The reaction dynamics of twisted intramolecular charge-transfer (TICT) state formation of 4-(N,Ndimethy1amino)benzoniitrile(DMABN) in methanol solvent has been studied theoretically. Molecular dynamics calculations were carried out for the SZstate of DMABN, and the reaction free energy surfaces were constructed as a function of solvation coordinate and torsional angle of the dimethylamino group. The solvation coordinate was defined by the potential energy difference between the S1 and S2 states. It was found that the potential barrier for the TICT state formation exists and the free energy along the solvation coordinate was well represented by the harmonic curves. To investigate the dynamics of reaction, the equations of motions for the solvation coordinate and the torsional angle were derived and stochastic trajectory calculations were carried out. The results were that (a) the reaction rate is significantly slower than the time scale of relaxation in the solvation coordinate, (b) the rate is about half of the transition-state theory prediction, and (c) the coupling between the motions along the torsional angle and the solvation coordinate is relatively small.
1. Introduction The formation of the charge-transfer state in solution is one of the fundamental processes in photochemistry. As is wellknown, the rate of this process is affected by static and dynamic properties of the solvent.' For example, in cases where the reaction is controlled by solvent motion, the rate constant is theoretically predicted to depend on the inverse of the longitudinal relaxation time of the solvent, Z L - ~ . Although some molecules reveal the charge-transfer rate depends on ZL-~,there is a series of molecules whose intramolecular charge-transfer rates have little correlation with the relaxation times of solvents.' 4-(NJV-Dimethylamino)benzonitrile(DMABN) belongs to the latter. DMABN exhibits dual fluorescence in polar solvents, and the formation of the twisted intramolecular charge-transfer (TICT) state in the excited state has been considered as the origin of the red shifts of the fluore~cence.~,~ The dynamics and energetics of the TICT state formation process of DMABN have been examined in many experimental works.4 In recent studies, Hicks et al. have proposed that the barrier height for the TICT state formation is controlled by the polarity of solvent in their experiments for the nitrile and alcohol ~ o l u t i o n s . They ~ ~ ~ have also pointed out that there is little dependence of the rate on the viscosity of solvent. Su and Simon have observed that the average lifetime of DMABN in the local excited state is shorter than Z L . ~ They have suggested that a fluctuation in the intramolecular coordinates is more important than the solvent diffusion. This conclusion was deduced from the analysis of their experiments on the basis of theoretical model developed by Sumi et al.839 To obtain a clear understanding of the mechanism of solutionphase chemical reactions, several statistical calculations such as Monte Carlo and molecular dynamics simulations for the solute-solvent systems have been performed by several researchers.10-20 In recent studies, Kat0 and Amatatsu carried out ab initio calculations on the excited states of DMABN and obtained pair potentials to describe the solute-solvent interact Present address: Department of Chemistry and Biochemistry, University of Colorado, Boulder, CO 80309-0215. Abstract published in Advance ACS Abstmcts, December 15, 1994. @
0022-3654/95/2099-0955$09.00/0
tions as well as the potential energy surfaces in the internal coordinates.20 The mechanism of TICT state formation in aqueous solution has been discussed on the basis of results obtained from Monte Carlo simulations. In the present work, molecular dynamics (MD) trajectory calculations have been carried out for DMABN in methanol solvent. We focus on the dynamics of TICT state formation in the S2 state of DMABN. In the TICT state formation process, the internal rotation of the dimethylamino group induces the charge transfer. We thus discuss the role of the solvent in the reaction on the basis of a microscopic reaction model by examining dynamics in both the internal rotation and a solvation coordinate which represent the fluctuation and reorganization in the solvent polarization. In the next section, the method of calculations and simulation results are presented. We first construct the two-dimensional reaction free energy surface as a function of the solvation coordinate and the torsional angle of the dimethylamino group. The solvation structures have been also discussed in terms of the radial distribution functions. In section 3, we introduce a reaction model to describe the dynamics of TICT state formation process. By deriving equations of motion for the torsional angle and the solvation coordinate, the correlation of the dynamics between them is discussed. To estimate the reaction rate, we have described the dynamics of the torsional coordinate with a generalized Langevin equation and stochastic trajectory calculations have been performed. We give a discussion on the mechanism of TICT state formation dynamics. Conclusions are summarized in section 4.
2. Molecular Dynamics Calculation 2.1. Method. We carried out the MD simulation calculations on a DMABN molecule in methanol solvent near the room temperature. We treated the solute and solvent molecules as rigid bodies, except for the two large-amplitude internal degrees of freedom of DMABN, the torsional angle z, and the wagging angle 8 of the dimethylamino group (see Figure 1). The potential energy of this system was evaluated with the parameters based on the ab initio calculations.20 The electronic 0 1995 American Chemical Society
Hayashi et al.
956 J. Phys. Chem., Vol. 99, No. 3, 1995
Y
41
$
Figure 1. Coordinate system of DMABN.
Y
v
structure of the S2 state in solution is described as a superposition of three diabatic states D,, Db, and D,. The D, state is of ionpair type and the Db and D, states are of covalent type?' The Hamiltonian matrix of the system is expressed as
H~~ = o,(e,t)
+ d,P;t(e,r&,r)
I)
I
G
"0
0
'I
2
O t
-12
(K,L = a , b , ~ (1) )
-6
-8
-10
Solvation Coordinate
-4
(kcal/mol)
s
where DKLis the energy matrix element between the K and Lth diabatic states and is the solute-solvent interaction energy for the state K . Energy of the Sz state in solution is calculated by solving the seqular equation
e
HC = ,IC
(2)
We approximated the solute-solvent interaction potential, which include the electrostatic and exchange-exclusionenergies, by a sum of pair potentials between the atoms of the solute and solvent. The methyl groups were regarded as extended atoms. The pair electrostatic potential was expressed as Coulomb interaction between effective point charges on the atoms. The effective charges on the atoms of DMABN were determined so as to reproduce the electrostatic potential of DMABN molecule. The pair exchange-exclusion potential which was computed by the Gordon-Kim mode121was fitted to the 12-6 Lennard-Jones function. For the solvent methanol, we used the parameters developed by Jorgensen.22 The details of calculations of these functions are given in ref 20, and the functions are summarized in Appendix A. When the MD simulation calculations are canied out, the Coriolis coupling between the internal vibrational motions and the rotational motion of the body-fixed frame of the DMABN molecule must be taken into account. We adopted the procedure given in ref 19 to include this effect. Integration of the equation of motion was performed with the five- and four-value Gear predictor corrector method, initiated by the Runge-Kutta method, for the first- and secondorder equation of motion, respectively. The number of molecules in the simulation box with the length of 29.4 8, was l DMABN and 216 methanol molecules. The mass density was set to be 0.759 g/cm3. We calculated the solvent-solvent longranged Coulomb interaction using the Ewald sum technique. We have implemented the potential tapering method for the computation of the solute-solvent interaction. The contribution of Coulomb interaction energy from the outside of the sphere was estimated with the reaction field approximation. Equilibrium MD calculations of 36 ps were performed, after 15 ps cooling and equilibration runs with a time step of 0.4 fs. Nonequilibrium MD calculations were also carried out with the same time step. No temperature control algorithm was used. The resultant temperatures are 298 f 14 K. 2.2. Free Energy Surface. MD calculations were performed to obtain the free energy surface on which the TICT state formation proceeds. We have chosen two coordinates to describe the free energy surface. One is the torsional angle of the dimethylamino group z, and the other is a solvation coordinate defined as follows. It is convenient to adopt the potential energy difference between neutral and ion pair states as the reaction coordinate in order to discuss charge-transfer reactions in polar s o l ~ e n t s . l ~We - ~ ~have therefore introduced
12
8
4
0
Solvation Coordinate
s
,
4,
!
20
10
16
(kcal/mol) I
30
Solvation Coordinate
s
I
40
(kcal/mol)
Figure 2. Free energy curves along the solvation coordinate at z = (a) 0", (b) 60°,and (c) 90". Symbols 0 and dashed curves are defined in text. The values o f t are indicated in each panels.
the solvation coordinate defined by the energy difference
s=
w,- w,
(3)
where W1 and W2 are the potential energy for the system with the SIand Sz state, respectively. It is noted that the SI state is of the covalent type20 and the S2 state is described as a superposition of the diabatic states D, D, as seen above. We calculated the free energy curves along the solvation coordinate s at seven points of fixed z spaced out 15" from 0" to 90". An additional calculation was carried out at z = 66". The free energy F(s) was computed with the probability density for the solvation coordinate P(s):
-
F(s) = -k,T In P(s)
+C
(4)
where C is a constant. The free energies at z = 0", 60", and 90" are presented in Figure 2, where the symbol 0 corresponds to the calculated values with eq 4. Assuming Gaussian distributions for s, the free energy curves can be expressed as
J. Phys. Chem., Vol. 99, No. 3, 1995 957
TICT State Formation of 4-(N,N-Dimethylamino)benzonitrile 2r
h
-
i 28
Y
v
,
I
I
1
,
I
1
1
I
,
I
1
1
O
x
?!
;-1 -21
0
'
'
,
30
60
90
Torsional Angle (deg)
Figure 3. Potential of mean force for the torsional coordinate. Symbols 0 are calculated values and a dashed curve is fitted one to the function Fc(t,4).
-20
0
k = @(s - s,n)2]-'
(6)
40
20
Solvation Coordinate
Figure 4. Free energy surface of the
s
S2
(kcal/mol)
state. Contour spacing is 1
kcal/mol. where p-' = ~ B T The . curves by eq 5 are also shown with the dashed lines in Figure 2, where we can see that the free energies are well represented by the harmonic potentials of eq 5. This indicates that the assumption of linear response to describe the motion along the solvation coordinate is reasonable at all the values of z. The potential of mean force for the angle z was also calculated with the importance sampling method and is shown in Figure 3. This potential was least-squares fitted to the analytic function Fc(zJV = 4): N
F,(T,N) =
C a ncos 2nz
c I
e
." * a
p
1.0
-
0.0
t
L * I n
.-Q
2.0
(8)
n=O
The profile of the potential curve along z in the solution is strongly modified from that in the gas phase where the energy of the Sz state monotonically increases along the torsional coordinate.20 The potential barrier for the TICT state formation appears in the present calculation. The top of barrier is located at z = 55" with a barrier height of about 1.8 kcal/mol. This is slightly underestimated compared to the experimental value of 2.6 kcal/mol measured in ethanol solution? The formation of TICT state is to be exothermic, AG = -1.7 kcalfmol, in the present calculation. The calculated value of AG is in good agreement with the experimental values of --1.5 kcal/mol measured in nitrile^.^ In this product well, the ion-pair-like diabatic state, D,, is the main configuration of the solute electronic wave function, as contrasted to the gas-phase case, because of electrostatic stabilization by solvation, whereas the D, state and the covalent type Db state are strongly mixed in the reactant well. At first, we estimate the reaction rate of the TICT state formation on the classical transition state theory (TST) with the one-dimensional curve in Figure 3 and the moment of inertia for the torsional motion defined in Appendix B. The resultant value of kTsT, 1.2 x lo-' ps-', is expected to be the upper bound of the actual reaction rate. The lifetime k ~ s ~ of - l 8.3 ps is too short compared with the survival time for the local excited state of 20 ps experimentally measured in cool methanol
solution^.^ The total free energy surface is given in Figure 4. It is of interest that the position of the minimum of curve in the solvation coordinate at a given value of t shifts to the larger value as t becomes larger in the barrier region (t = 45-66'). It is noted that the force constant for the solvation coordinate decreases with increasing z.
I
4.0
Distance
6.0
8.0
(A)
Figure 5. Radial distribution functions between the centers of mass of DMABN and methanol at t = O", 6 6 O , and 90".
Before proceeding to the TICT formation dynamics, we discuss the solvation structures by examining with the solutesolvent radial distribution functions (RDFs). The RDFs between the centers of mass of DMABN and methanol at z = 0", 66", and 90" are shown in Figure 5. We can observe a large peak centered at about 5.3 A for all the values of T. At z = 90" a new sharp peak at 3.6 8, appears, indicating that there exists strong coordination of CH30H to the N1 and C1 atoms of DMABN for the TICT state. This strong coordination is clearly demonstrated by the N1-0 and C1-H RDFs given in Figure 6a,b, respectively. We can also observe a strong hydrogen bonding between the terminal N2 atom of DMABN and the H atom of CH30H as in Figure 7 where the narrow peak at 2.0 A is found. 3. Reaction Dynamics
3.1. Reaction Model. In the previous section we calculated the free energy surface as a function of two coordinates, the torsional angle t and the solvation coordinate s. There is, however, a very large number of degrees of freedom for the solvent motion which affects on the dynamics of the TICT state formation. We have dealt with the complicated motion of the solute and solvent molecules by projecting them onto z and s. As seen above, the linear response relation seems to be satisfied for the motion along s at all the values of z. We therefore attempted to use the harmonic bath model for describing the motion along the solvation coordinate. The potential energies of the SIand S 2 states are expressed as
958 J. Phys. Chem., Vol. 99, No. 3, 1995
Hayashi et al. from the origin at x* is 5=x-x
E
." a n ." Y
*
(11)
Taking into account the fact that the X* depends on z, the Hamiltonian of the Sz state is written as
1.0
r)
6 l/$?b22 where
0.0 2.0
4.0
6.0
Distance
I
8.0
g(z) = g,
(A)
+ ozx*
(13)
'
cI
0.0
+ g'(t)x + ws(z) (12)
2.0
4.0
6.0
J
Si(Z)
(A) Figure 6. Radial distribution functions between (a) the N1 atom of DMABN and the 0 atom of methanol and (b) the C1 atom of DMABN and the H atom of methanol OH part at t = On, 66",and 90".
=dS
-
Distance
zi=O
li....;
zj=0
(15)
2.c
E
."0 Y
where S represents the length measured along the direction
s =8 4
(16)
As is easily seen, the solvation coordinate eq 3 is proportional to the coordinate S defined in eq 16:
0.0
' ''
..
2.0
4.0
6.0
(A) Figure 7. Radial distribution functions between the Nz atom of DMABN and the H atom of methanol OH part at t = 0", 66",and Distance
s = p-1'2s (17) where p is regarded as the effective mass for the solvation coordinate because p satisfies the equipartition theorem 1
90".
W&z)
+
= 1/&f32x g&>x
+ 9At)
(K= 1 , 2 ) (9)
where x represents the fluctuation in the potential energy of the system, i.e., the bath coordinates. As was previously shown,lg the solvation coordinate, eq 3, can be related to the displacement vector measured from the minimum energy crossing point between the two surface. The minimum energy crossing point of the two surfaces x* is given by 1 1 n i(z) = X
*
kBT --
(18)
where the dot denotes the derivative with respect to time. Figure 8 shows the effective mass p evaluated MD calculations with fixed z. The effective mass is monotonically decreases along the torsional coordinate. Because the effective mass depends on z, we have used the mass-weighted solvation coordinate S = p*I2s, instead of s to construct the Hamiltonian. This provides us with more simple expressions of the equations of motion. The bath coordinates are partitioned into the massweighted solvation coordinate and the other X whose directions are orthogonal to S, i.e.
i
2= &
(19)
pi = 0
(20)
where
(10) for each value of z and the displacement of bath coordinate 2
By applying the canonical transformation from the old variables
J. Phys. Chem., Vol. 99, No. 3, 1995 959
TICT State Formation of 4-(N,N-Dimethylamino)benzonitrile I
x
4w
0.01
4
,
I
30
60
90
I
o
"
"
"
"
'
I
Torsional Angle (deg)
Figure 8. Change of the effective mass ,u along the torsional coordinate. Dots are calculated values and a dashed curve is fitted one to the function F,(t,3). (Z#) and (zp,) to the new ones (S,Ps), (Xp),and (Z,Pr), the
Hamiltonian is rewritten by
1
0
,
,
,
,
, , > 1000
,
'
I 2000
Frequency (cm-1)
T = '/?P
+
1 'I2P: 4--{Pr 21
- P,,
+
- x*"(DP SPs)}2
Figure 9. Components of the reaction coordinate vector Si at t = 45" and 66". The values of t are indicated in each panel.
the velocity autocorrelation function: (SS(t)) = ( k B r / p ) x s ;cos(0,t)
(33)
i
where Pcorrepresents a Coriolis term:
P,,, = S S W
+ ?PSPs
(24)
In eq 23,SZ is regarded as the frequency of the harmonic free energy curve along the solvation coordinate: Q
=
fi2
- ii, PQ
2aQQ-22 OQP
where -2
- -1
2-
0 -SOS
thp; = aQ;= up$ WPQ,2= (0; - Q2)Zi
(aQ;)jj= (U0Qp2Qrj= fipp,2dij WQQ,2 = 0;6,
+ (fi2 - 0;- O;)z,Zj
and Vr is potential energy which only depends on z:
vr = -11 252-2-2g (z> + $(@ where
g = (g'O-2g)-'[g;O-2S
e)]
- p(@ -
The present MD calculations with fixed z provide several quantities in the Hamiltonian eq 21. Components of the direction vector S are computed with the Fourier transform of
We found that the components Si were very similar for all the values of z. The results at z = 45"and 66" are shown in Figure 9. In these figures, the components Si can be characterized by four different bands. We can observe two distinct bands B1 and B2 centered at 350 and 700 cm-', respectively. Two broad bands B3 and B4 are also seen in both low-frequency, 100200 cm-', and high-frequency, 800-1200 cm-', regions. The B2 band is originated by the librational mode of methanol OH g r ~ u p ' and ~,B ~ 1~band by the vibrational motion of the wagging angle 8, respectively. We can also roughly estimate the origin of the other bands with the calculations of the Fourier components of the velocity correlation functions for various valiables. In Figure loa, the components of kinetic energy of the solute molecule at t = 45" are presented. These are two bands corresponding to the B1 and B3 bands. The B3 band is therefore the contribution from the translational and rotational motions of the solute molecule. Figure 10b shows the components of the total kinetic energy of the solvent molecules, which we can see has a broad band in the high-frequency region corresponding to the B4 band. We calculated the frequency SZ with eq 25, and the results are summarized in Table 1. The values of 52 are almost independent of z. We also calculated the frequencies with the alternative relation
Q,=Ji&
(34) The frequencies Q1 are also included in Table 1. It is found that the frequency calculated by two different methods are very similar to each other. Hereafter, we use the average value of SZ1, 185 cm-', as the frequency. In deriving the equation of motion from the Hamiltonian eq 21, we assumed that the components Si and 0~are set to be constants because the distribution of components Si has a similar form for all the values of z. The Coriolis term P,,, in eq 22 thus vanishes. The Hamiltonian with this approximation provides us with a set of equations of motion for the coordinates X, S, and z. It is straightforward to obtain the equations of
Hayashi et al.
960 J. Phys. Chem., Vol. 99, No. 3, 1995
' A '
'
'
'
'
'
'
1 ..... 01
E
rlG(t)
2.0
i
g
1.0
.C
Y
E
0.0
I 0.0
0.1
0.4
0.3
0.2
0.5
Time (ps)
Figure 11. Friction kernels ~ ( tat) t = 45" and 66" and sc(t). 1.0
Y
- ETCF r = 45' ___ ETCF r = 6G0
I
0
1000
OLE
2000
Frequency (cm-1) Figure 10. Fourier components of a velocity correlation function of the kinetic energy of (a) the solute molecule and (b) the solvent molecules at t = 45".
TABLE 1: Freauencies Q and QI t (deg) P (cm-l) PI (cm-l) t (deg) P (cm-') 0 207 179 60 191 15 170 184 75 205 30 209 197 90 203 45 187 180
.........,.
9 1
(cm-') 188 187 179
motion for z and S by projecting out the bath variables 2. If z is fixed, we can obtain an equation which represents the dynamics along the solvation coordinate:
(35) Note that eq 35 is the so-called generalized Langevin equation (GLE). Vs is the harmonic free energy along the mass-weighted solvation coordinate:
v, = '/,S2*(S - s,n)2 The friction kernel
v ( t ) is
(36)
0.0 0.0
0.2
0.4
.....1 . - - ..__..... ,
0.6
0.8
1.0
Time (ps)
Figure 12. Comparison of relaxation in the solvation coordinate between the equilibrium time correlation function (ETCF, eq 43) at t = 45" and 66" and the prediction by overdamped Legevin equation (OLE,eq 42).
in the solvation coordinate as seen below. Accordingly, we can approximately describe the relaxation dynamics with the Langevin equation in which the friction constant is calculated by i j = JvG(t) dt = (&/2)avG(O)
(41)
where we approximated the friction kernel by vG(t). Moreover, it seems to be valid that the inertia term of the Langevin equation is neglected because q/2Q is larger than 1. After all, the solvation dynamics is described by the over-damped Langevin equation:
expressed as
+
dV , ijs = - - F,(t) dS and satisfies the fluctuation-dissipation theorem with the random force Fs(t):
(F,Fs(t)) = v(t)k,T
(39)
The friction kernels at t = 45" and 66" computed with eq 38, are displayed in Figure 11. These friction kernels decline very rapidly and are followed by small oscillations with frequencies mainly corresponding to the B2 band in Figure 9. The initial decay curves are well fitted to a Gaussian function:
which were also drawn in Figure 11. The time constant a was almost the same for all the values of z and was calculated to be 19 fs. This time scale is much shorter than that of the relaxation
In Figure 12, the relaxation in the solvation coordinate predicted by eq 42 is presented. It shows exponentially decaying behavior with time constant qlQ2 of 0.15 ps. We have compared the prediction of eq 42 with the normalized time correlation function of the solvation coordinate obtained by the fixed z, 45" and 66", equilibrium MD calculations:
C(t)= (SS(0) ss(t))/(ss(o)2)
(43)
where 6s is the fluctuation in the solvation coordinate, 6s = S - (3. The Correlation functions are also shown in Figure 12. Although we have introduced several crude approximations, the time character of the relaxation predicted by eq 42 is in reasonable agreement with that of the time correlation functions in the short-time region.
J. Phys. Chem., Vol. 99, No. 3, 1995 961
TICT State Formation of 4-(N,N-Dimethylamino)benzonitrile If the GLE eq 35 is replaced by the over-damped Langevin equation eq 42, the equation of motion for t is represented by the following form:
~10-5
- g',fOY( t ) -
wo
__--_ _ - -_ _ - - -
.._.--+ __..-
0.0 0
+
+
, 60
30 T
(degree)
Figure 13. Comparison of the dispersions of 9%(0)and Fcalculated by eq 67 (+) and eq 48 with the approximation eq 69 (a dashed curve) at t = 0 66". A solid curve is fitted one to the function Fc(t,3)for
,..
(am.
w, = -1/gt2"-2g2
+@
(47) and R(t) and Y(t) are given in Appendix C. The right-hand side of eq 44 consists of five terms. It is easy to show that WO is the potential energy along the torsional coordinate and coincides with the free energy curve presented in Figure 3. The f i s t and second terms may be regarded as the frictional and random force terms of GLE, respectively, since g can be approximated by a linear function of t as seen below. The random force fit) satisfies the fluctuation-dissipation theorem:
.....
0.0
0.2
0.4
7
= 4Y
T
= 66"
EXP
0.6
0.8
1.0
Time (ps)
The time constant of the friction kernel, ijlQ2, coincides with that of the relaxation in the solvation coordinate. This means that the relaxation and fluctuation in the solvation coordinate described with eq 42 is coupled with the motion of t as the frictional and random forces. The third term represents Coriolis forces due to the shift of the origin of the bath coordinate x* as the function of t (see Appendix C). Y(t)in the fourth term is given by
Y(t) = &r)
+ thQQ-2thQ;S(t)
(49)
and satisfies the relation
(50)
(SSY(t)>= 0
The fourth term therefore represents the forces from the bath modes perpendicular to S. It is shown in eq 44 that the dynamics of t is influenced by the various terms through gz and @ - @. Unfortunately, we could not deduce these variables from the present MD calculations. We have therefore described the motion along t by following GLE:
A = -W, - Z
sf,((?-
t') t(t') dt'
+ %(t)
(51)
and directly obtained the random force s(t)from the MD calculations with the fixed t. The friction kernel [ ( t ) is estimated with the fluctuation-dissipation theorem:
(%q(?)> = zg(t)kBT
Figure 14. Normalized auto correlation functions of 5at 0" and 60". An exponential (EXP) decay curve is fitted one.
regions. The normalized autocorrelation function of .3$at t = 0" and 60" are presented in Figure 14. The time behavior of the friction kernel is almost invariant in the product and barrier regions. These correlation functions decay rapidly and are well fitted to a Gaussian function with time constant of 64 fs. 3.2. Stochastic Trajectory Calculation. We have attempted to estimate the reaction rate of the TICT state formation with stochastic trajectory calculationsz5 based on eq 51. It is, however, difficult to apply the stochastic trajectory methods to the GLE directly because this requires very large matrix inversion procedures in order to compute the random force which satisfies the fluctuation-dissipation theorem. Thus, we have approximately split eq 51 into a couple of equations:26
zi:= -y'(z)y - W r
sz2
1
5 - -5y ( t )
Y =-y
1 + -F,(t) 5
(53) (54)
where the second equation is an over-damped Langevin equation and the relation (Ffy(f)>
= 2EkBTS(t)
is satisfied. As seen in section 3.1, if y ( t ) is linear to equations are equivalent to eq 51 with
(55) t,these
(52)
(ss,
Figures 13 shows the dispersion of 3, at t = 0-66". The dispersion increases as T becomes larger and hence the frictional and random forces become stronger with approach to the barrier region. In the product region t = 75-90", the dispersion is much larger than that in the reactant and barrier
(57)
Hayashi et al.
962 J. Phys. Chem., Vol. 99, No. 3, 1995 I
bI
2ot 10
I
'
a&
i 20
0
40
80
60
Time (ps)
Figure 16. Distribution of the average survival time in reactant region. Dots are calculated values and dashed curves are fitted one to an
exponential function.
,
1.0
-
7
, = 45"
6 is determined by
fitting an exponential decay curve to the autocorrelation function of 4,and this experimental curve is presented in Figure 14. y ( t ) is determined by the relation
At first, we calculated the relaxation in the torsional coordinate with both the stochastic trajectory and nonequilibrium MD calculations. About 30 initial configurations were sampled with at least a 1 ps interval from the fixed-z calculations. We computed the root-mean-square displacement and the time profiles are shown in Figure 15. They are characterized by a Gaussian initial decay with the time constant of about 0.15 ps. The results by the stochastic trajectory are consistent with those by the nonequilibrium MD calculations, although we introduced various approximations to calculate the dynamics of z with the stochastic trajectory calculations. We employed the stochastic trajectory calculations in order to estimate the reaction rate. For the distribution of the initial state, we used the equilibrium distribution of the harmonic potential:
I
6,
Kqv= '/2k,z? k, = I(n$AG*~ST)2 and about 10000 configurations were sampled with this potential. Figure 16 shows a distribution of the survival time of each trajectories. This represents the time profile of the population in the reactant well and is well fitted to the exponential curve, which is also presented in Figure 16, with the time constant t,Qof 16 ps. This average survival time is larger than ~ T S T -and ~ the transmission coefficient defined by = (tRRICTST)-l
(62) is 0.53. We also computed the Grote-Hynes transmission coefficient K G H : ~ ' = zJob
I
0.2
0.4
0.6
0.8
1.0
Time ( p s )
Figure 17. Normalized time correlation functions ( 6 S S ( t ) ) @ S a at t = 45",66", and exponential decay (EXP) calculated with eq 46.
is about 0.5, and hence the frictional and random forces acting on z by the solvent motion are important in describing the reaction dynamics. In the present treatment, we have separated the solvent motions into two contributions; one is the motion along the solvation coordinate, and the other is the bath motion projected out in deriving the GLE (eq 35). To examine the importance of each contribution, we have analyzed the forceforce correlation function (Ss(t)). The dispersion of the correlation function is already shown in Figure 13. As seen in eq 46, the direct force from the solvation coordinate motion is given by
(3s
where k, was estimated from kTST:
KGH
0.0
L ,
(63)
where wb is the barrier frequency and &z,) is the Laplace transform of ((2) at the barrier top. The calculated KGH, 0.59, is in good agreement with the present stochastic trajectory result. 3.3. Implication of Reaction Dynamics. The stochastic trajectory calculations revealed that the transmission coefficient
9 f t ) = j'dS(t)
(65)
and it is straightforward to show that the correlation function has the form
where
(m
=(6sa2Q2/kBT In deriving eq 66, we used the relations ( S S q ( t ) ) = (SS9ft)) = (kBT/Q2)j'e-'n2'7)t
(67)
(68)
We have computed the cross-correlation function (dSs(t)) by the MD calculations, and the results are given in Figure 17. It is shown that the initial decay of (dSs(t)) is well represented by the exponentially decay curve with the time constant of +jlQ2 = 0.15 ps, which supports the relation eq 68. The dispersion of Fobtained with eq 67 is included in Figure 13. Compared with (3s is Figure 13, the direct contribution of the solvation coordinate is much smaller, which indicates the random force acting on z is mainly originated from the bath motions.
J. Phys. Chem., Vol. 99, No. 3, 1995 963
TICT State Formation of 4-(N,N-Dimethylamino)benzonitrile Moreover we can alternatively estimate introducing the approximation for g: g' = -Q2stw-2sol-1/2)' + Q2st0-2g'1
~
(mwith eq 48 by
supported by the Grants in Aid for Scientific Research from the Ministry of Education.
-Qzs'o-2sol-l/2>'
Appendix A (69) where SZ2Sfw-2Zis computed to be 0.9. It is noted that the 2 calculated with eq 69 is almost linear to z in the barrier region. The resultant dispersion of F i s presented by the dashed line in Figure 13, which also shows that the direct contribution of S is quite small. These analyses tell us that the direct coupling of fluctuation of the solvation coordinate with the dynamics of the torsional motion is small and the motions of these two degrees of freedom are mainly driven by thermal fluctuation of the bath degrees of freedom.
4. Conclusion In present work, we have carried out the MD trajectory calculations on DMABN in methanol solvent to elucidate the mechanism of TICT state formation. The two-dimensional reaction free energy surface was constructed as the function of the torsional angle z and the solvation coordinate S. As for the solvation coordinate, we have chosen the difference of potential energies between the S1 and S 2 states. The free energy curves along the solvation coordinate are well represented by the harmonic curves. It was found that the barrier height of the potential of mean force for the torsional angle is about 1.8 kcal/ mol. In the barrier and product regions, the D, state, which has the ion-pair-like electronic structure, becomes the main configuration of the solute molecule. It is therefore expected that the barrier height depends on polarity of solvent as suggested by Hicks et al. in their experimental work^.^^^ To investigate the dynamics of the TICT state formation process on the free energy surface, the Hamiltonian of the system was constructed by applying the harmonic bath model to describing the motion along the solvation coordinate, and we derived the equations of motion for the torsional angle and the solvation coordinate. The motion in the solvation coordinate is well described with the overdamped Langevin equation (eq 42) and the relaxation time was calculated to be 0.15 ps. The relaxation and fluctuation in the solvation coordinate are coupled with the equation of motion for the torsional angle as the frictional and random forces. We reduced the dynamics of z to a GLE, and the stochastic trajectory calculations were performed in order to estimate the reaction rate of the TICT state formation. The resultant survival time of 16 ps is in good agreement with the experimental value.' This survival time is much longer than the relaxation time in the solvation coordinate, and hence the reaction is little controlled by the dielectric relaxation time of the solvent. We also evaluated the transmission coefficient for the reaction. The resultant transmission coefficient of about 0.5 represents a marked deviation from the transition-state theory prediction. It was found that there is little direct coupling of the random force acting on z with the fluctuation in the solvation coordinate by computing the correlation functions between them. Thus the fluctuation of the bath degree of freedom seems to play an important role in the reaction dynamics.
Acknowledgment. The authors are grateful to Professor T. Okada and Mr. K. Nishiyama for their valuable comments and discussions. Numerical calculations were carried out at IMS Computer Center and the Supercomputer Laboratory, Institute for Chemical Research, Kyoto University. This work was
The potential energy surfaces of the S1 and three diabatic states were described by the analytic function 1
wAt,8) = C(a,,
+ an202+ an606)cos 2nt
(70)
n=O
The off-diagonal matrix elements of the D1 state with the Dz and D3 state were represented by the function
D,J(t,8) = (b, cos z
+ b3 cos 3z) cos 8
(71)
and the D2-D3 off-diagonal element was by D2,(t,8) = (bo
+ b2
COS 22
+ b,
COS
4t)(l - c - c COS 8) (72)
respectively. The DMABN-methanol intermolecular potentials were given by the sum of Coulombic and exchange-exclusion terms. For the Coulombic part, the z and 8 dependence of the effective charge on the atoms in DMABN was represented by the function
+
qi(z,8)= (ai,o ai,2cos 22
+ ai,4cos 42 + ai,6cos 62)( 1 +cos
e) (73)
where the subscript i denotes the atoms or extended atoms of DMABN. The values of parameters in these functions are available from the authors.
Appendix B When the motion of z is discussed, the moment of inertia I must be defined. The moment of inertia is expected to depend on a distribution of the wagging angle 8. To estimate I, we have defined the center of inertia for z and 8 at first. If DMABN is located in the reference coordinate x, with the N1 atom at the origin, the N1-CL bond along the z axis and the aromatic ring on the x-z plane. The center of inertia coordinate xc is defined by x, = ux,
(74)
U = CBA
(75) (76)
0
c=
[rSisin,,s)
(77)
R'
sinz,
COST
where 8 R and t~ represent the centers of inertia for 8 and z, respectively, defined by (79)
964 J. Phys. Chem., Vol. 99,No. 3, 1995
where and are moments of inertia of the benzonitrile moiety and two methyl group for the motion of 8, and frBN and ILM are the analogous for the motion of z. Kinetic energies for 8 and z is expressed as
P=
d”:)
G=
where pe and p r are the conjugate momentums of 8 and z, respectively, and mi denotes the mass of the atom i. p r is derived from eq 81 as follows: (84) where gy are elements of G-l. By taking an average over the distribution of 8, the second term in right-hand side of eq 84 vanishes. The moment of inertia I is hence given by (85) I = (
