J. Phys. Chem. 1083, 87, 2881-2889
H-’=
[
]
(ef)x
(egIx ( e h I x
(ef)s
(eg)s ( e h I y
(ef)Z
(eg)s
(A-17)
(eh)r
Appendix B. Definition of the Projection Factor P
The projection factor p , as defined in eq 8, permits a decision concerning the sign of no,by correlating the direction of translation of helix 2 in the (fg)plane, given by t, = t,ef tgeg (Figure 2c) and the direction of the axis f ’ (Figure 2a) around which the axis of helix 2 is rotated with respect to the axis of helix 1 by the angle /3 (Figure 2b). According to the definition of the angles a,/3,y in the Euler rotational operation, the axis f ’lies in the (fg)plane (Figure 2a). The angle between axes f and f ’is a. Looking along f’ toward the positive direction of the f’ axis, /3 is defined to be positive if the axis of helix 2 is rotated in the clockwise direction relative to helix 1 (Figure 2b). In order to define the sign of Q,, we consider it as a dihedral angle, formed by three vectors, eH and eHp(eq 3 and 4), pointing in the positive direction dong the axes of helices 1 and 2, respectively, and D which is perpendicular to both eH, and eH2(defined here to be directed from HI to Hz).According to the usual definition of the sign of dihedral angles,12this sign is taken to be positive if the unit vector that is farther from the observer is ro-
+
2881
tated clockwise with respect to the closer unit vector, or conversely, if the nearer one of the two unit vectors is rotated counterclockwise with respect to the farther one. The sign is negative for opposite rotations. In order to relate !doto 0, it is necessary to decide which helix is farther for an observer looking along the f ’axis in the positive direction. For this purpose, we consider the direction of the translational displacement t, relative to the direction of the f ‘axis. If e, is a unit vector along the f ’axis and e, is a unit vector pointing in the direction of the translational vector t,, as defined in eq 6 and 7, respectively, then the projection of e, on e, is given by eq 8. If p is positive, then the absolute value of the angle between f’and t, is 90° and helix 2 is moved toward the observer, so that it becomes the nearer helix. In this case, /3 and nohave opposite signs. The relationships described result in eq 9. If p = 0, t, is perpendicular to t,, i.e., the two helices are at the same distance from the observer, measured along the f ’axis. Rotation by /3 around f’ takes place in the plane defined by eH and t,, i.e., the two axes of the two helices are coplanar. fn this case, the sign of Oo cannot be defined uniquely. Registry No. Ac-(L-Ala)lo-NHMe, 85251-49-6.
Time-Correlation Functions from Computer Simulations of Polymers Thomas A. Weber’ and Eugene HeHand Bell Laboratories, Murray Hi//, New Jersey 07974 (Received: December 14, 1982; I n Final Form: March 2, 1983)
Many of the rapid relaxation processes in polymers are related to conformational transitions of bonds from one state to another. Earlier research, based on Brownian motion computer simulation and kinetic theory, revealed that conformationaltransitions occur as single events and also as cooperative pairs. That research focused on kinetic rate parameters, but the experiments generally measure time-correlation functions of the system’s properties. This paper is a report on the results of computer simulations of polymer time-correlation functions. The decay with time of both initial conformational states and initial vector orientations has been studied. The correlation functions fit well to a two-factor functional form, suggested by recent theoretical research. The first factor is an exponential decay arising from single transitions. The second is a diffusional term related to the correlated pair transitions. The total transition rate when compared with the earlier kinetic studies is in reasonable agreement. Since the functional form has been shown here to fit the simulations, it is suggested that it will be of value in interpreting relaxation experiments.
I. Introduction Fundamental to the rapid relaxation processes which occur in polymers are the conformational transitions from one rotational isomeric state to another. NMR, dielectric relaxation, dynamical light scattering, and fluorescence depolarization are a few of the experimental techniques which may be used to study such relaxation processes. The maximum amount of information which can be extracted 0022-3654/83/2087-288 1$01.50/0
from these experiments would be a time-correlation function of some orientational quantity (generally a vector or a second-rank tensor). This orientation might be defined, for instance, by an electric dipole, a polarizability anisotropy, or the direction of a proton relative to a 13C nucleus. Because the direction is fairly rigidly defined relative to the polymer backbone, and reorients when bonds undergo conformational transitions, one could hope 0 1983 American Chemical Society
2882
The Journal of Physical Chemistry, Vol. 87,No. 15, 1983
to unravel the mechanism of conformational transitions from such measurements. In practice this has proven to be difficult. It appears to be more promising to proceed in the other direction; i.e., hypothesize a mechanism and deduce from this a variety of experimental implications. This a priori approach is useful in two ways: (1)it suggests functional forms to which data may be fit to extract relevant parameters; and (2) when a variety of data fits the theoretical predictions, then this may be regarded as an experimental verification of the assumed mechanism. In this paper we assume the role of experimentalists; i.e., we take the results of computer Brownian dynamics simulation of polymer motions and generate time-correlation functions. We determine not only vector correlation functions but also conformational-state time-correlation functions (defined in section 111). Although the latter are experimentally inaccessible, they are more firmly linked to theory and help to bridge the gap between the proposed mechanisms and orientational time-correlation functions. The underlying mechanism which we assume for conformational transitions is one that emerges from prior studies-the simulations by Helfand, Wasserman, and Weber,l and the kinetic theory analysis by Skolnick and Helfand.2 The fundamental question in a consideration of conformational transitions in polymeric backbones is how rotations can take place about a bond without being overwhelmingly resisted by the need to swing the polymeric tails attached to that bond. The answer is that, as torsion occurs in the transforming bond, in particular as the bond's torsional angle traverses the transition barrier, there is an accompanying distortion of neighboring degrees of freedom. These distortions are such that asymptotically the tails remain stationary. One type of distortion is particularly important in localizing the motion. If the first neighbor bond of the transforming bond is trans, then the second neighbor is parallel. Cranklike counterrotation of this second neighbor bond is then very effective for localizing the motion. After the torsional angle of the transforming bond is well over the energy barrier, the distortions, on average, must dissipate. Frequently, however, one finds that the second neighbor bond is brought all the way to the point of a transition of its own. Thus, in a sequence of conformational transitions, one observes both individual transitions
t t g" (trans s gauche plus and trans cranklike pair transitions ttg" 2 g'tt
(1.1) 2
gauche minus) and (1.2)
ttt 2 g*tg' (1.3) Hall and Helfand3 have derived the form of the conformational time-correlation functions for systems undergoing similar transitional processes. It is on the basis of an approximate application of these arguments that the Brownian dynamics simulations reported here will be analyzed. The present research is based on Brownian dynamics simulations of polymers with trajectories run 10 times longer than our earlier studies.l We find that this does not improve in any major way the rates that we previously reported based on hazard analysis. These longer runs, however, were essential for the calculation of meaningful correlation functions. (1) E. Helfand, Z. R. Wasserman, and T. A. Weber, Macromolecules, 13, 526 (1980). (2) J. Skolnick and E. Helfand, J . Chen. Phys., 72, 5489 (1980). (3) C. K. Hall and E. Helfand, J . Chem. Phys., 77, 3275 (1982).
Weber and Helfand
TABLE I: System Parameters moperts
value
P = S-lm
1.0x 105
m
0.014 0.153 2.5 X l o 9 70.53 1.3 X 10' 1.0 1.3108 -1.4135 -0.3358 2.8271 -3.3885 6.634 X l o 5 12.36 2.933 44.833 - 60 -120
bo 7bIm 0 0
YO lm
00
a, a2 a3 a4
a,
vglm
E* E, Ecis
o*
Og
'
units ns- ' kg mol-' nm J ns- kg- ' deg J kg-'
J kg-' kJ mol-' kJ mol" kJ mol-' deg deg
A preliminary report of a few of the results of this paper has been presented earliera4 In the previous paper only the data for the two higher temperatures were available. In addition, an improved model for the cooperative transitions is included here. The plan of the paper is as follows. In section I1 the polymer model is briefly described and the method used to generate the trajectories is discussed. In section I11 the definition of the various vector autocorrelation functions and the trans and gauche conformational time-correlation functions are given. In the next section models are introduced to determine relaxation rates for the processes of eq 1.1-1.3 from the correlation functions. The relaxation functions proposed by Hall and Helfand3 are found to fit both autocorrelation and cross-correlation simulation data exceedingly well. A comparison is made between the transition rates determined from the time-correlation functions and rates determined more directly by hazard analysis (a statistical study) of times between transitions. The vector autocorrelation functions are considered in section V. 11. Molecular Model
The interaction potential used to model the n-alkanes and infinite molecular weight polyethylene has been described in detail elsewhere .l However, a brief review will be presented here for completeness. The polyethylene skeletal backbone is represented by united-atom carbon centers, Le., hydrogen atoms collapsed into the carbons. This carbon backbone undergoes motions which can be described as bond stretchings, angle bendings, and torsional rotations about the bonds. The carbon-carbon bond vibrational potential for bond i (the bond between vertices i - 1 and i) of length bi,is taken as Vb(bi) 1/yb(bi - bo)' (2.1) where Yb is the force constant and bo is the bond length at 0 K. (The parameters used in this polyethylene simulation are listed in Table I, and correspond to potential A of ref 1.) The angle Bi a t vertex i between bonds bi and bi+l is kept near the tetrahedral value by the potential
vo = '/Zys(COS
Bi - COS 80)'
(2.2)
where cos Bo = The potential used in the Ryckaert and Bellemans5 n-butane simulation is employed to govern the torsional motion: (4) T. A. Weber, E. Helfand, and Z. R. Wasserman, in press, (5) J.-P. Ryckaert and A. Bellemans, Chem. Phys. Lett., 30,123 (1975).
The Journal of Physical Chemidty, Vol. 87, No. 15, 1983
Time-Correlation Functions
v,
5
c ak k=O
= ,y@
COSk
4i
(2.3)
The trans configuration is chosen as the zero of energy for this potential. The trans-gauche barrier is E* = 12.36 kJ/mol; the energy of the gauche minima is Eg = 2.933 kJ/mol; and the cis barrier is Ecis= 44.833 kJ/mol. No other carbon-carbon interactions have been explicitly included in this model. Most notable excluded volume interactions, such as the pentane effect which virtually prohibits g*g‘ pairs, are neglected. Also carbons widely separated along the chain may pass through one another (a phantom chain model). We do note believe that the presence of such spurious motions affects in any important qualitative way the mechanics on the short time scales considered here. The infinite molecular weight polyethylene chain is constructed by periodic replication of a segment consisting of 200 units. This additional boundary condition eliminates end effects and makes each carbon center equivalent. The position of a carbon center 200 units removed from any vertex, i, is ri+200
=
ri
+
(P201
- r1l)ez
(2.4)
where e, is a unit vector in the x direction and ( lr201 - rll) is the average end-to-end distance of a 200-unit chain. Equation A.8 of the Appendix gives an excellent approximation for ( lrzol- rll ) for the potential specified by eq 2.1-2.3. When the initial coordinates of all 200 carbon centers and the temperature are specified, the Langevin equations of motion for high-friction Brownian particles p(dr,/dt) = -(l/m)ViV
+ Ai(t)
(2.5)
are solved to generate a representative trajectory. In eq 2.5, m is the mass of the carbon center, /3 is the friction constant divided by the mass, and V is the total interaction potential composed of the bond length, angle, and torsional parts. The incoherent effect of the medium (solvent or other polymer molecules) surrounding the carbon center is replaced by the Gaussianly distributed random forces A,(t) which have zero mean and covariance (Ai(t) A,(t?) = (2pk,T/m)GijS(t - t ? l
(2.6)
In eq 2.5 and 2.6 one sees that p and t enter only in the combination tip. Thus, changing friction constants amounts only to rescaling time; i.e., all rates are proportional to l/p. Numerical trajectories were generated by using Helfand’s extensiong for stochastic differential equations of the Runge-Kutta method and working to second order in the time step. The integration time step was At = 5 X lo4 ns. A t each temperature in the range 330-425 K the system was equilibrated by running a trajectory until a few thousand transitions had occurred. This trajectory was then continued for from 2 x IO6 to 5 x lo6time steps until 40 000 transitions had occurred. A typical calculation required on the order of 10 h on the CRAY-1 computer. 111. Time-Correlation Functions
In a previous series of calculations’ the transition rate as a function of temperature was determined for the Brownian dynamics polyethylene simulations using hazard analysis.’ It was not possible to extract decay rates from
2003
autocorrelation functions generated from these trajectories because the total length of time of the simulation was too short to reduce noise. By increasing the length of the trajectories by a fador of 10 it is now possible to determine transition rates from the autocorrelation functions. In this study, we report on several vector and bond conformational autocorrelation and cross-correlation functions at three different temperatures. In a previous paper4we have given a preliminary report on the bond and one of the vector autocorrelation functions at two temperatures. Vector autocorrelation functions are defined as N
( h ( ~ ) . h ( O ) ) = (1/N)
J h i ( ~ ’ + ~ ) - h i ( ~dT’/JdT’ ’) i=l
(3.1)
where hi(7) is some unit vector associated with either vertex i or bond i at time T . The summation runs over all 200 bonds or vertices and the integral is from 0 to a sufficiently large time (in practice the length of the simulation less the length of time over which the correlation function is calculated). Four different orientational autocorrelation functions, as defined below, have been calculated: a bond autoconelation function associated with bond i; and chord, out-of-plane, and bisector autocorrelation functions associated with vertex i. The bond autocorrelation function is calculated from the unit vector along bond i bi = bi/bi (3.2) The chord autocorrelation function measures the loss of correlation of a vector which we have found to be a good indicator of the direction of the polymer backbone. The chord vector is the unit vector connecting the midpoints of two adjacent bonds, i.e.
The out-of-plane vector associated with vertex i is defined as (3.4)
i.e., it is the unit vector perpendicular to the plane of bi and bi+l. Finally, a vector in that plane and nearly perpendicular to the chord vector is defined by the bisector vector of the angle between bonds i and i + 1 (3.5)
The bond, chord, out-of-plane,and bisector autocorrelation functions at 330 K are shown in Figures 1-4, respectively. Notice that the bond and chord autocorrelation functions do not reach equilibrium due to the slowness of relaxation of long wavelength vectors in polymers (see section V). Combined plots of the four different vector autocorrelation functions are shown in Figures 5-7 for the three temperatures studied. Besides knowing the vertex coordinates at each time, we also know the bond’s conformational state, i.e., whether a bond is trans (t),gauche plus (g+), or gauche minus (g-1. Conformational autocorrelation and cross-correlation functions have been constructed as follows. Define the function As,i(T)
(6) E. Helfand, Bell. Syst. Tech. J.,58,2289 (1979). (7)N. R.Mann, R. E. Schafer, and N. D. Singpurwalla, ‘Methods for Statistical Analysis of Reliability and Life Data”, Wiley, New York, 1974).
= 1 if bond i is in conformational states at time = 0 otherwise
7
(3.6)
2884
The Journal of Physical Chemistry, Vol. 87, No. 15, 1983
Weber and Helfand 10
08
0 6
r 3
0 4
32
got -3 2
00
o 4
8
6
-0 2
time i n s )
time (ns) Figure i. Bond vector autocorrelation function at 330 K defined by eq 3.1 and 3.2. The dlamonds represent the simulation data. The solid curve is the best fit to the data using eq 5.1.
Figure 3. Outaf-plane vector autocorrelation function at 330 K defined by eq 3.1 and 3.4. The diamonds represent the simulation data. The solid curve is the best fit to the data using eq 5.1.
I 4
2
0
2
6
4
6
time (ns)
t i m e in:) Flgure 2. Chord vector autocorrelation function at 330 K defined by eq 3.1 aml 3.3. The diamonds represent the simulation data. The solid curve is the best fit to the data using eq 5.1.
Flgure 4. Bisector vector autocorrelation function at 330 K defined by eq 3.1 and 3.5. The diamonds represent the simulation data. The solid curve is the best fit to the data using eq 5.1.
Henceforth, we will use the shorthand s = -, 0, + for g-, t, g+, respectively. The conformational time-correlation functions to be studied are defined by
in our study. Furthermore, we will use the known symmetries to reduce the statistical uncertainties c,(l,7) = C s ( - l , 7 ) (3.9) C + ( 1 , 7 ) = c-(1,7) (3.10) Le., we report 1/2[c8(l,7)+ cB(-l,7)] and c*(1,7) = 1/2[c+(l,7) + c - ( l , ~ ) ]Finally, . note that ( A s , L ( 0is) )just the average fraction of bonds in state s PsE (AS,LO)) (3.11)
c , ( ~ , T )=
(As,r+l(7)As,L(0)) -
([As,b(0)12)-
(As,L(0))2 (3.7)
The averaging is performed as a time average over the trajectory, and, since all bonds are equivalent, we also average over i : 1
N-1
(A8,r+d7) A 8 , r ( 0 )=) - C N t , r=o
tM
d7’ As,L+1(7’+7) As,L(7’) (3.8)
with t M taken as the length of the run less the largest
7
and since
takes only the values 0 and 1 ( [As,,(O)J2)= Ps
(3.12)
The cross-correlation functions, such as co&(1,7),may be calculated by using the normalization, symmetry, and time reversal properties of the g* and t states once the corre-
The Journal of Physical Chemistry, Vol. 87, No. 15, 1983
Time-Correlation Functions 1
1
1
1
2885
out-of-plane c
2 1
I
1-;oo.o9 1
01 00
10
05
I
0
15
time (ns)
00
0 5
10
15
time (ns)
Figure 5. Vector autocorrelation functions at 330 K. The diamonds represent the simulation data. The soli curves are the best fits to the data using eq 5.1.
Flgure 7. Vector autocorrelation functions at 425 K.
c
a 0
o u t -of - p l a n e
T
= 330
K
T
= 372
K
z1 la d
0 30
10
05
15
time (ns) Figure 6. Vector autocorrelation functions at 373 K.
lation functions defined by eq 3.7 are determined. The trans autocorrelation function cO(0,7), and the correlation function for second-neighbor bonds, c0(2,7), have been calculated for the three different temperatures of the simulations and are shown in Figures 8 and 10, respectively. The gauche auto and second-neighbor correlation functions, c+(o,7)and ct(2,7), have also been calculated and are shown in Figures 9 and 11, respectively. We have previously p ~ s t u l a t e d lthat - ~ the most important mechanisms for conformational transitions are
ttt
e k'if
g*tg'
(3.13~)
1
2
time ( n s ) Flgure 8. Trans conformationalautocorrelation function, c &O.T), defined by eq 3.7. The diamonds represent the simulation data. The solid curves are the best fir to the data using eq 4.4 and the constraint on Xl/Xf1 (see Table IV).
with rates as indicated adjacent to the arrows. The forward and reverse rates are related by detailed balance. Reaction 3.13a represents the rate for single, independent transitions. The other two reactions (3.13b and 3.13~)are the most important in which cooperative transitions are involved as can be seen in Table 11. Two transitions are considered to be cooperative if they occur within time steps of one another. Specifically any two transitions which occur within 5 ps of one another for the 372 K simulation are considered as cooperative. The cutoffs for the other two temperatures are scaled by a Boltzmann factor. When one uses this criterion, over 29% of all transitions are found to be cooperative (a slight overestimate since statistically some pairs of independent ones are expected within T ~ " ~ ) .Over 80% of all cooperative transitions are of the type specified by reactions 3.13b and 3.13~.The overall rate at which transitions occur, A, is
Weber and Helfand
The Journal of Physical Chemistry, Vol. 87, No. 15, 1983
2006
1 T
I
330 K
=
0 5 1 L
i ci
-
‘r
-
530 h
T
=
372 k
T
=
325
i.
P i Y
3 “i
w
0.0
K
C
1 I
0
J 1
0
2
2
time ( n s )
time (ns) Figure 9. Gauche conformational autocorrelation function, c*(O,r), defined by eq 3.7. The diamonds represent the simulation data. The solid curves are the best fit to the data using eq 4.4 and the constraint on h,Ih‘, (see Table IV).
Flgure 10. Trans conformational cross-correlation function, c 0 ( 2 , r ) , defined by eq 3.7. The diamonds represent the simulation data. The solid curve is eq 4.7 using the parameters of Table I V which fit the autocorrelation function data.
given by a sum over all allowed reactions of the product of the number of bonds transforming in the process, the rate of that reaction, and the equilibrium probability of the initial state: X = 2kofPo + 2ko,.P* + 4k’lfPo3 + 4P$*[k’1Q* + ( k ~ f+ k1r)PoI (3.14) (where P , P- = P+). Let us define for the three processes, cy
A, =
(3.15) 2 P , = 1 we can write
(kafkar)lj2
Using detail balance and Po +
T
=
372 K
T
=
425 K
.c M
kor =
k’1f
[
>]1’ 12),o-
(3.16b)
Po
1 - Po =2P0 A’ 1
(3.16~)
2P0 ktIr = -A’, (3.16d) 1 - Po which implies that the overall transition rate is h = 2Po(l - PO)[A& + 2PO(h1 + A’I)] (3.17) where 2
=
0000.
112
[ Po(l - P o ) ]
(3.18)
One method by which the overall transition rates can be measured is hazard analysis.’ The probability that a bond undergoes a transition between time T and r + d r is h ( r )d r where h ( r )is the hazard rate. Throughout this analysis, transitions have been determined by a first passage time criterion. The three regions, g-, t, g+ are delimited by the maxima of the torsional potential. A transition is regarded as having been completed when the barrier is crossed and the torsional angle reaches the minima of the potential well of another region. The transition time is taken as the time between transitions.
0
I
time (ns) Flgure 11. Gauche conformational cross-correlation function, ch(2,r), defined by eq 3.7. The diamonds represent the simulation data. The solid curve is eq 4.7 using the parameters of Table I V which fit the autocorrelation function data.
The hazard is defined as H ( r ) = Jrh(r’) d7’
(3.19)
and the (large T) slope of this function is the transition rate. Figure 12 shows the hazard plots constructed from the first passage time data for the three different temperatures. (The method for preparing hazard plots from data is discussed in ref 1.) The total transition rates determined in this way are reported in Table IV for the three temperatures. This will be compared with rates extracted from the correlation functions. It should be mentioned that the rates estimated from the slope of the hazard plot are something of an overes-
The Journal of Physical Chemistry, Vol. 87,No. 15, 1983
Time-Correlation Functions
correlation functions for this model may be determined exactly to be3
TABLE 11: Cooperative Transitionsa temp, K % cooperative T cut
g'tt $ ttg' t t t "t g' tgF t t t 2 g'tg' g'tt 2 ttg' gFg't "t tg'g'
tg't
fgfgtg'
gTg't
"t tg'gF
g'g't Ztg' tg'g' f g'g't tg't "t gTg'g'
"1
330
372
425
29.3 1649 5513 4151 438 411 3 94 352
29.2 5394 4162 3 89 344 428 376
29.4 607 5188 4163 330 303 569 442
192
24 8
313
179
24 1
331
76
79
90
40
41
39
1000
c+(Lt) E 4(-4+,i+l(t)A+,i(O)) - 1 = exp(-2Xot) exp(-2Xlt)Il(2Xlt)
+
a Two transitions which occur within considered to be cooperative.
rCUt time
30
steps are
/
,/' =
425K,,/'
20
2 21
6
c
10
00
0 1
0 2
time ( n s ) Figure 12. Hazard functions, see eq 3.19, for the three temperatures studied.
timate. This is due to correlations between transitions, which manifests itself also in a slight curvature in the hazard plots. We plan to analyze this in greater detail in the future. IV. Conformational Correlation Functions From a theoretical point of view it is easier to understand configurational correlation functions than orientational correlation functions. Therefore, we begin with an investigation of the former. Hall and Helfand3 have discussed the theory. They approach the problem by way of solutions of a master equation. Before studying the polymer model considered in the present paper, they consider a simpler model which has two, rather than three states per bond. These bond states, which they label p j = fl (i = 1, ..., N), are considered as equally likely at equilibrium. The allowed kinetic processes are (i) single bond transitions between states + * (4.1) at the rate X, and (ii) pair transitions in which one bond's index is changed in an opposite way from its nearest neighbor's +- * (4.2) at the rate XI (this is the analogue of counterrotation). The
-+
(4.3)
where Il is a modified Bessel function. The result is easily understood in physical terms. The initial state of the system corresponds to knowledge of the state of bond i, the state of all other bonds being random. This bond is returned to randomness by the single processes of eq 4.1, which accounts for the factor exp(-2Xot). It can be shown that the pair process of eq 4.2 transfers the undecayed knowledge to a neighbor (i - 1or i + 1) and leaves the state of i random. Therefore, the probability of the state of bond i 1 being is enhanced above the random value by the probability that a one-dimensional lattice random walker, taking steps to the right or left at a rate A,, is displaced exactly 1 positions a t time t. This is given by exp(-2Xlt)ll(2Xlt). Hall and Helfand next consider a three-state problem with processes such as eq 3.13. (However, they do not require that the middle bond be trans. Also they take all states as having equal energy. We will return to these points.) This problem no longer reduces to a Markoffian random walk problem. The reason is that there are more ways to exchange with a neighbor which is trans (which can go up or down) than with one of the gauches. Thus, as the knowledge of the original state transfers, knowledge of the preferred state of the other site is retained at site i. A simple approximation adopted by Hall and Helfand is to drop this correlated information. This leads to correlation functions containing modified Bessel functions, as in the simpler case. In their paper Hall and Helfand only consider explicitly the case of three states of equal energy, although they outline the method for unequal energy. They have since performed the calculation for a gauche state having an energy E , greater than trans.* Also, to account for the fact that the bond separating second nearest neighbors must be trans for counterrotational transitions, eq 3.13b or 3.13c, to occur X1 and will be multiplied by Po and thus the explicit state of the center bond need not be considered. This is again in the spirit of ignoring multibond correlations. The result is
+
/
T
2887
+
Co(117)
=
exp{-[XoE + 2PdX1 +
X'1)17)I(1/2,1[2~0'0(~'l- x1)7I
(4.4)
for 1 an even integer, and no correlation between the state of bonds which are odd-integer neighbors. One can understand the appearance of two Bessel function terms in ct as arising from the different relaxation rates for the two modes of deviation of the three states from equilibrium. One is a symmetric mode in which the trans is favored over its equilibrium population at the expense of the two gauches equally (or vice versa depending on the sign of the term). The other mode is an antisymmetric one in which one gauche is favored totally at the expense of the other. Only the former mode can contribute to co in the present approximate scheme. (8) E. Helfand and C. K. Hall, unpublished result.
Weber and Helfand
The Journal of Physical Chemistry, Vol. 87, No. 15, 1983
2888
TABLE 111: Four-Parameter Fit to Eq 3.7 temo. .,K
330
372
4 25
Po,fit A,. ns-'
0.6091 0.4771 4.2328 0.5217 0.1232 0.7530 5.1865
0.6037 0.6976 8.2672 0.7925 0.0959 0.7716 9.6347
0.5585 1.6433 12.7275 1.1426 0.0898 0.8024 15.9880
,: ns
1
P o ~ ' ns-' I, / A ])fit (h',/h,),,,,~ation A, ns @'I
TABLE V: Vector Autocorrelation Function (Acf) Fitting Parameters 313 425 330 temp,K Po,simulation
KO K l
a li
TABLE IV: Three-Parameter Fit to Eq 3.7 temp, K
330
372
425
0.6114 0.4368 2.3947 1.8031 4.5916 5.8 .
0.6072 0.6029 4.4012 3.3959 8.2716 10.3
0.5606 1.5405 6.6346 5.3238 13.9454 17.1
-
P0,fit
ns-' &A., ns P,h ,, ns A , ns h h a z a d , ns \o,
KO
K 1
01
li
K 1
0.5586
0.7715 8.6133 0.1741 0.1415
Bond Acf" 1.1236 13.3690 0.1694 0.1263
1.8646 20.7900 0.1678 0.3932
0.6617 5.4201 0.2718 0.1537
Chord Acf" 0.9336 7.9051 0.2555 0.1299
1.5780 12.2249 0.2566 0.4051
Bisector Acf" 3.0157 4.5331 14.0252 23.3100
KO
7.2046 35.2113
Out-of-Plane Acfa
To begin with, the trans and gauche autocorrelation function data (1 = 0) from the simulations were simultaneously fitted to the functional forms, eqs 4.4 and 4.5. The corresponding to the best values of Po, Xo, X1, and functional fit are reported in Table I11 for each of the three temperatures. (It should be pointed out that Po is determined primarily by its explicit appearance in the equation ( A , , L [ ~ ( A,,L[~L(0)I) ~)l = Ps(l - Ps)c,(o,7) + P:
0.6006
0.6110
(4.6)
rather than from the appearance of Po in eq 4.4 or 4.5.) These parameters do not yield a good fit to the secondneighbor correlations c0(2,7) and ck(2,7). Also it is found that the four parameters can be varied considerably without seriously degrading the goodness of fit. It thus appears that four parameters give too much flexibility to the fit; i.e., the uncertainties of the parameters are so great that the best-fit values have little significance. In particular, it is difficult to pin down the ratio X1/Xfl, which is very important for the 1 = 2 correlation function fits. Therefore, it was decided to take this raio directly from the simulations. It is easy to show, by using detailed balance, that the ratio of the geometric mean of the rates of reaction 3.13b to that of reaction 3.13c, X,/X',, is equal to the ratio of the number of occurrences of these reactions. These are compiled in Table 11. Again one sees that the parameters of Table I11 do not fit the ratios. Good fits were then made to the cO(O,s) and cf(0,7) with the three free parameters Po, &, and XI, and with A', specified to match the observed ratio. The solid curves of Figures 8 and 9 show the fits to the data with the parameters of Table IV. Moreover, the same parameters now fit the 1 = 2 correlation data fairly well, as is shown in Figures 10 and 11. In particular, these parameters give a very small trans correlation, c0(2,7), in excellent agreement with the data, whereas the parameters of Table I11 give much larger values. It is significant that, although the agreement to the autocorrelation functions was obtained by fitting, yet the second-neighbor correlation curves are also in excellent agreement with the data. The overall rates of conformational transitions are also listed in Tables I11 and IV. These are to be compared with rates determined from the long-time portion of the hazard plots shown in Figure 12. Maximum likelihood, straight-line fits were made to the hazard plots. Half the data points, those corresponding to the shortest times, were censored so as to eliminate the curvature of the hazard plots a t small time. As was mentioned earlier the slight long-time curvature cannot at present be eliminated so
3.0321 11.8765 0.0018
KLl K1 CY
a
The units of
K
~ k,
4.0917 21.1864 0.0006 and
i(
7.5415 30.8642 0.0000
are ns-'.
there is an uncertainty, perhaps 2070, in the results. The rates determined in this way are listed in Table IV. The agreement between rates determined by hazard analysis and from correlation functions is fairly good.
V. Orientational Correlation Functions Our theoretical knowledge about the functional form of the orientational correlation functions is much slimmer. Preliminary work indicates that there are many more factors to be taken into account and the results, in detail, will be much more complicated, unless the many factors combine to yield in some approximate sense a simpler result. At this stage we will adopt this hopeful approach and attempt a functional fit to the simplest function which incorporates the elements that we have found to be most important. In particular, we shall assume that there are individual reorientation processes which lead to an exponential loss of orientational correlation at one rate. Also there are cooperative processes which contribute to a loss of correlation modified by the probability of the undegraded information diffusing 1 positions in time 7. This leads to a factor e - * ~ ~ l ~and ( ~ ~again 7 ) , we will assume only one such process. However, this is not enough, because we know that certain molecular vectors (the bonds and chords) are correlated with long wavelength vectors like the end-to-end distance, which decays very slowly. Thus, we will assume that a part of the vector correlation is slowly decaying, and for simplicity we use a single exponential relaxation. In total then we will fit orientational autocorrelations to the functional form @(7)
= (1 - 0')eXp(-Ko7) eXp(-K17)lo(K17)
+
CY
eXp(-p7) (5.1)
The use of a complementary error function, Le., replacement of e x p ( - ~ ~ r ) Zby ~ (e2*1' ~ ~ erfc ~ ) [ ( 2 ~ ~ 7 ) "has ~ ] ,been suggested from theoretical considerations by Valeur, Monnerie, and J a r r ~ .The ~ two functions agree a t large time but the erfc form could not be made to fit the important curvature of the autocorrelation functions when K ~ was T O(O.l) as well as the loform. (We have not attempted fits to the functional forms suggested by Jones (9)B. Valeur, L. Monnerie, and J.-P. Jarry, J. Polym. Sci., Polym. Phys. Ed., 13, 675 (1975).
The Journal of Physical Chemistry, Vol. 87, No. 15, 1983
Time-Correlation Functions
and Stockmayerlo or Bendler and Yaris.") The values of K ~ K, ~ a, , and p obtained from a nonlinear least-squares fit of the simulation data to eq 5.1 are given in Table V. The solid curves of Figures 5-7 show how exceedingly well this functional form fits the data, perhaps better than the conformational correlation function fits. As expected, the parameter a is quite small for correlation of vectors perpendicular to the chain (when a is small, p cannot be determined). One caveat which must be needed when attempting to extract rates from the fits is that, since the initial decay goes like exp[-(no + K ~ ) Tit] ,is somewhat difficult to determine the individual parameters. I t is our belief that a better functional form than eq 5.1 can be found to fit the orientational correlation function. Work in this direction is in progress.
Appendix Average End-to-End Distance. The average end-to-end distance between two carbon centers separated by n bonds is related to the mean-squared distance by ( Irn - rol)= (8/3r)lI2((rn- ro)2)1/2 (A.l)
sf3
2889
0
(-4.4)
where cg
= (cos 0)
Sg
= (sin 0)
c4 = (cos 4 )
(A.5)
Equation A.3 may be further approximated for the large-n case as
(Ir, - roI2) = N b 2 ) + 2 N b ) 2 [ ( T ) / ( 1- (T))Il1 (-4.6) For a bond constrained by the potential of eq 2.1 (A.7a) (A.7b) where (6 = kBT/Ybb:. distance squared is
Therefore, the average end-to-end
where r, - ro is Gaussianly distributed. The vector connecting these two points is N
rn - ro = C bi i=l
(-4.2)
where bi is the ith bond vector. Therefore, for a potential which is a sum of independent bond length, bond angle, and torsional angle terms the mean-squared end-to-end distance is
Cg
(Ir, - rOl2) = C(bibj) ij
n-1
= n(b2) + 2n C ( b ) ( T ) k l l ( b ) k=l
The integrals of eq A.5 are calculated numerically for the potential functions eq 2.2 and 2.3:
(-4.3)
where ( b ) is an average bond length and (T) is an average rotation matrix which rotates bond i into bond i + 1. The average rotation matrix is (10)A. A.Jones and W. H. Stockmayer, J.Polym. Sci., Polym. Phys. Ed., 16,847 (1977). (11)J. T.Bendler and R. Yaris, Macromolecules, 11, 650 (1978).
=
JrCOS
0 Sin 8 eXp[-Vg/k~T) dO/
