Faraday Discuss. Chem. SOC., 1987, 83, 199-211 Brownian Dynamics of Chain Polymers Marshall Fixman Department of Chemistry, Colorado State University, Fort Collins, Colorado 80523, U.S.A. The effect of constrained bond lengths and angles on the Brownian dynamics of wormlike polymer chains is calculated analytically and compared with dynamical simulations. An analytical correspondence between constraints and internal friction is developed and illustrated. The approximate numerical evaluation of the constraint matrix and its correlations, characteristic of previous work, is replaced by an equilibrium simulation. Good agreement is found between the constraint formalism and dynamical simulations. With fluctuating hydrodynamic interaction included, the constraints generate a coupling between local motions and translational diffusion which causes a small decrease in the translational diffusion constant.The effects are similar to those found in a previous phenomenological study of the effects of internal friction on translational diffusion and do not vanish with increasing chain length. It is well known that the classical Gaussian model192 of polymer chain dynamics is inadequate for short, stiff chains or for local motions that affect the moderate and high-frequency properties of flexible chains. The necessary refinements in the model seem fairly obvious for the systems considered here: dilute solutions at the 0 point, Instead of the very flexible bond lengths and bond angles allowed in the Gaussian model, only modest vibrations should be permitted and motion in the dihedral angle subspace should be restricted by rotational barriers.Progress with these improvements has been slow in each of the two approaches that promise a wide generality of application. On the one hand, a fundamental approach involving the representation of the diffusion operator on a complete basis ~ e t , ~ ' ~ although successful for very short gets bogged down in numerical difficulties for chains with more than a few bonds. An obvious analytic approximation involves the pre-average of an inverse constraint If bond angles as well as lengths are constrained, this approximation fails qualitatively for the local motions of long chains. On the other hand, the same improve- ments in predictive power which are the goal of fundamental approaches are easily accomplished by the addition of internal friction to Gaussian models." The manoeuvre can be formally justified as an approximation to memory integrals in Zwanzig-Mori projection but the actual applications depart considerably from this solid foundation.An attempt is made here to supply some of the missing structure.I3 Our confidence that a useful link can be developed between formal theory and internal friction goes back to Brownian simulations on chains with constrained bond lengths and angles and barriers to internal r0tati0n.l~ Virtually all of the normal mode structure of Rouse-Zimm theory seemed to be preserved by the more realistic model, and the similarity to the predictions of internal friction models was noted. But the results were sketchy and the connection was not made in a systematic way.A more extensive treatment has now been started. The objective is to determine first whether the analytic treatment of constraints is successful and practical for the treatment of freely rotating chains with arbitrary stiffness and length, initially without rotational barriers, and whether there is a natural correspondence to the treatment of such chains by the addition of internal friction to Gaussian models. It will be seen that these questions can be answered affirmatively. 199200 Brownian Dynamics of Chain Polymers To some extent the technical difficulties are circumvented rather than confronted. The required averages and correlations involving the constraint matrix are obtained from equilibrium simulations rather than by clever approximations.In retrospect, however, the averages are seen to have simple properties; they converge rapidly with increasing chain length to intensive parameters characteristic of the backbone structure and friction constants. Basic Equations A very brief summary of the constraint formalism will be given. Earlier work4 has been modifiedI3 by the addition of internal friction forces to the basic Langevin equation. The purpose of this addition is to allow a translation of constraints into the language of internal friction. At no point in the actual calculations is internal friction included as an independent contribution. The Langevin equation with neglect of inertia has the symbolic form m d V / d t = F H + F 1 + F P + F R --+ 0 ( 1 ) where FH is the frictional drag force due to the relative motion of beads and solvent, F' is the frictional drag due to the relative motion of beads with respect to each other, FP is a potential force and FR is a random force.These symbols all stand for column matrices of order 3( N + 1). The bead positions R,, i = 0, . . . , N, and velocities V, are likewise written as 3( N + 1 ) dimensional arrays R and V, respectively. where 1I' is the distribution function of R. In the simulation algorithm" which will be referred to below, FR is a stochastic variable. The frictional drag force is F H = -Po( V- v"), where Vo is the unperturbed velocity of the solvent and Po is a friction matrix. Po= (Pi1 + T ) - ' , where Pd is a diagonal array of bead friction constants and T is the matrix of hydrodynamic interaction HI (e.g.the Oseen interaction, although in practice the Rotne-Prager interactionI6 has always been used because it keeps Po positive definite for all configurations). The internal frictional force F' = - Qo V and the coefficient matrix Qo will be taken independent of R in order to make simple contact with memory integrals. The constraining forces have known direction (determined by the form of the backbone vibrational potential) and form one part of the total potential force F'. The unknown magnitudes are determined by the criterion that bond lengths and bond angles remain constant in time. The details of the calculation may be consulted The result is that In analytical work FR= -kTV In v = G ( p , V + F " + F R ) (2) where Fs is the remainder of FP after constraint forces are subtracted, and G= G-GMAG, G=(p,+Qo)-' and MA=AMAT. The matrix M is identical to that constructed b e f ~ r e .~ MAG is a projection matrix that prevents the bead motion from violating the constraints; the inverse of M, not M itself, is known explicitly. A,,, i, j = 0, . . . , N, stands for the elements of a constant matrix. It arises from the transforma- tion to bond vector coordinates b = ATR ( b , = R, - except for bo which is the centre of mass); AT is the transpose of A. G is a diffusion matrix in eqn (2). In the absence of constraints or internal friction the diffusion matrix would be H = P i ' . It is seen that constraints and internal friction affect the diffusion matrix in a superficially different way.If only one of these effects is present then the former leads to a subtraction from H, and the latter leads to an addition to H - ' . The subtraction is a major one if bond lengths and angles are both constrained, and approximations to M can lead to unphysical results such as negative eigenvalues in the diffusion matrix. The effect of internal friction is saved from that danger by the easily satisfied requirement that Qo is non-negative, i.e. that it has non-negative eigenvalues.M. Fixman 20 1 Introduction of the drift velocity into the conservation equation for 9 gives a diffusion equation for ZIr. This equation is usefully written in terms of @ = */ZIre, where qe is the equilibrium distribution function: With a scalar product defined as ( g , , g2) = (glg2),, i.e.as an equilibrium average, 2 is self-adjoint and for any functions g1 and g,. The explicit form of the source term S ( t ) is known in terms of V(' and external f o ~ c e s , ~ . * ~ but will not be needed here. The development of quasi-analytical approximations for transport coefficients is generally based on the representation of 3, S ( t ) , and @ ( t ) on a truncated basis set. Alternatively, if the transport coefficient has already been expressed in terms of time correlation functions of dynamical variables yk( t ) , y k ( t ) = exp ( - 2 t ) y k , then it suffices to represent 3 and the yk on the truncated basis set. We will refer here to the Zwanzig-Mori projection formalism3.12 rather than the Fixman-Kovac truncated rep- resentation~~ as a theoretical basis for this application, because the former provides an explicit remainder in terms of memory integrals.However, the memory integrals will be discarded after a brief contrast between the ways in which memory integrals and internal friction correct the truncated formulae, and the formalisms then become equivalent. The column array of y k ( t ) is designated y ( t ) . The generalized Langevin equation is X? d y ( t )/ d t = - LT, y ( t ) + K ( T ) r o l y ( t - T ) d T + A( t ) ( 5 ) 0 where I'o = ( yyT)e and eqn (4) gives Lij (yiTYj)e= kT((VTYi)G(Vyi)>e- (6) The kernel K ( T ) is a time correlation function of random forces A which are not correlated with y. Without the memory terms eqn ( 5 ) is a truncated matrix representation of the diffusion equation, as used in Fixman-K~vac.~ Eqn ( 5 ) is applicable to the caluclation of a correlation matrix I?( T) = ( y ( ~ ) y ~ ( 0 ) ) , .Improvements on dT( t)/dt =: -LI';'y( t ) + A( t ) or r( T ) = exp (- Lr&-)ro ( 7 ) which are obtained from eqn ( 5 ) after suppression of the memory integral, mist be based on an expansion of the set of yk or on an empirical treatment of the memory kernel. The imposition of constraints is one e ~ c e p t i o n , ' ~ and amounts to an evaluation of the time integral of K ( T ) over a short time interval that encompasses the vibrational relaxation times and a restriction of the yk to those variables which are essentially constant over such an interval. The use in eqn (6) of the constrained diffusion matrix G rather than the unconstrained G implies that the part of K ( T) representing the decay of vibrational forces has been assumed to have infinitely rapid decay and can be written as S ( T ) & , where KO is a constant matrix.These forces then modify the dissipative factor L in eqn (7), which is replaced by L - K O , rather than the equilibrium factor r i l , and in this respect behave similarly to internal friction. The first question to be addressed is the adequacy of eqn (7) in its application to a standard worm model of polymer chains with a minimum set of y k .202 Brownian Dynamics of Chain Polymers The Y k will be chosen as Rouse modes or as the direct products of such modes. The Rouse vector modes are defined as N qi = C Q.-b. j = O ?I J where QOo = 1, QOj = Qjo = 0 for j > 0, and An additional transformation b = ArR allows the q values to be expressed in terms of bead coordinates.For vector-mode correlations the basis set is composed of Y k = q f , 1 C k S N, where qz is the x component of qk. The tensor basis set consists of the products yk = q;qi. Time correlations of such modes suffice to describe dielectric relaxation, if the dipoles are rigidly affixed to the backbone, and stress relaxation, if HI is pre-a~eraged.~”~ However, that sufficiency does not guarantee than eqn (7) is accurate, i.e. that the basis set is large enough that memory integrals are negligible. This adequacy is judged by a comparison with dynamical simulations. The matrix elements of To are just ( ~ ~ y , ) ~ and their evaluation from equilibrium or dynamical simulations is clear in principle, once the y1 and the model are defined.For the calculation of L,J from eqn (6) it is desirable to transform from the original column matrix V y J to gradients with respect to Rouse modes, since the yJ have a simple dependence on them. The transformation to b derivatives, V = AVb, gives L I J = kT((v ’YI lTB(’ bYJ i > e (10) where B = ATGA = B - BMB; B = ATGA. (11) Except for the possible inclusion of internal friction Qo in G = (Po+ QO)-’ and the retention of the centre of mass coordinate bo, B and M are as defined previ~usly.~ For the orthogonal transformation from b to q coordinates, we use V h = QVq, where Q is the matrix of Qkl and V f = d / d q , . With B4 = QBQ, L I J = kT((V4r,)TBq(VqyJ)),’ (12) Approximations and Internal Friction We briefly summarize the major simplifications of the constraint formalism that have been tested by a comparison with results from the full formalism or with dynamical simulations.Several conclusions are anticipated. ( a ) Suppression of Constraints The suppression of constraints means that the constraint matrix M is discarded; i.e. G is replaced by G, or B is replaced by B in the expression for L, eqn (12). Suppression of constraints and internal friction restores the theory to the Bixon-Zwanzig formalism in a version that neglects memory integrals; this is very nearly a Gaussian formalism. The results are qualitatively reasonable for mode correlations but not for stress relaxatim. The results are always quantitatively poor for short, stiff chains, and are worse for short-wavelength, rapidly relaxing modes than for long-wavelength modes.The failure for rapidly relaxing modes does not diminish with increasing chain length. However, the low-frequency motion of long, flexible chains does not require retention of con- straints.M. Fixman 203 ( b ) Diagonal L Approximation In this approximation the off-diagonal elements of the matrix representation of 3 on the Rouse basis set are discarded. This will be seen to be an excellent approximation under almost all circumstances for uniform worm models. It is worse f,, non-uniform chains. When the approximation is valid the formalism and the results ma, ' - e interpreted in terms of internal friction, as will be shown below.We expect that this approximation will remain valid for chains with rotational barriers. (c) Factored Tensor L Approximation A related approximate factorization of the diagonal elements of L formed on the basis of tensor modes will be seen to be reasonably accurate. For the tensor modes eqn (12) We define the factored approximation to Lkk as A superscript v on L or Lik 3 kT( B;Ik), and coupled to ( q i ) 2 , then Lkk = Lfkk will be a good approximation. = 2 k ~ ( B 3 , ( ( ~ ; ) ~ ) , = x;,r;,, . (14) designates an evaluation on the vector basis set. Thus = (( qi)2)e. If fluctuations in the constraint matrix are weakly ( d ) Internal Friction Description This is not presented as another approximation. Our purpose is to show that the constrained diffusion matrix without internal friction can be replaced by the uncon- strained B with internal friction, and that internal friction has simple properties in Rouse coordinates if approximations ( b ) and ( c ) are valid. The internal friction matrix Qo appears in combination with the ordinary friction matrix Po= H-' in the matrix elements L,j of the time evolution operator. Specifically the diffusion matrix B in bond vector coordinates is, from its definition, where Bo= ATHA and A.is the 3 ( N + 1) x 3 ( N + 1 ) matrix A o = A-'QoAT-'. If Qo generates a frictional force on bead i which is due to beads i* 1 and which is equal to A ( V l i l - V ; ) , as in Bazua-Williams,'* then A. is A times a unit matrix in any representa- tion. More generally, eqn (15) gives where A:= QA,Q and B: = QBoQ We will assume that A:, although not generally proportional to the unit matrix, is to a sufficient approximation diagonal with elements A&.This is also a well known approximation' that we use for B,9; the diagonal elements are designated B&. Consequently, eqn (16) associates a friction constant with each Rouse mode. It remains to be demonstrated that internal friction of this form can emulate constraints, i. e. that BY with internal frictjon can approximate B4. With the substitution of BY from eqn (16) for BY in eqn (12) for L, and the use of a pre-averaged, diagonal B,9 and a diagonal A; in eqn (16), the factored approximation is always valid and we may write, in analogy to eqn (14), B = ATGA = AT( H - ' + &)'-'A = (B;' + Ao)-l (15) Bq QBQ = (Bg-' + A:)-' (16) Superscript I on L indicates use of the internal friction approximation, and superscript v on L or To indicates the vector as against the tensor basis set.It is now evident that204 Brownian Dynamics of Chain Polymers 0.0 n 0 h v 0 -1.0 > 0 % v U E - - 2.0 100 2 200 Fig. 1. Comparison of simulation and analytical tensor mode correlation functions for a 10 bead chain with b = 10, constrained bond angles equal to 160°, d s = 5, mode k = 5 , and no HI. L,/p = 0.6. ( a ) Dynamical simulation, ( b ) diagonal L approximation, (c) full analytic result, ( d ) omission of constraint matrices, ( e ) constraints and off -diagonal elements in the equilibrium correlation matrix To are omitted. the h z k can be chosen to make L& = L L k or L L k = L k k , and that both equalities can be satisfied simultaneously if the factored approximation is adequate.In summary, the adequacy of the diagonal L approximation guarantees that the effects of constraints can be imitated by an internal friction matrix which is diagonal in Rouse coordinates. The adequacy of the factored L approximation guarantees further that these diagonal elements are the same for vector and tensor relaxation properties. Results We begin with a summary of conventions. All of the graphical and tabular results are expressed in reduced units. The unit of length is a reference bond length brer, the unit of energy if kT, the unit of time is P,,,bf,,/kT, where Pref is a bead friction constant given by 3 rrqObref. Friction constants are specified through the equivalent Stokes diameter d as ratios to bref.No special notation is used to designate quantities expressed in these reduced units. Chain stiffness is measured by the ratio LJp, where L, is the contour length and p is the persistence length. Results obtained from the truncated matrix representation, eqn (7), will be described here as ‘analytical’ in order to distinguish them from the results of dynamical simulations. Comparison of Results with Dynamical Simulations Dynamical simulations were performed only for freely rotating chains without HI. The algorithm used for constrained simulations, with or without HI, was recently de~cribed.’~ Results are shown only for a nine-bond chain with L J p ==: 0.6. Quite similar comparisons were found for five-bond chains with larger and smaller values of LJp.Fig. 1 showsM. Fixman 205 Table 1. Full and approximate analytical results for a constrained bond angle chain with an angle of 163", b = 5, and d s = 3 (pre-averaged HI is included, L,/p =: 2); N = 40 ( a ) tensor mode friction constants ~ k P k 1 0.258 x lo3 0.229 x lo3 0.292 x lo2 0.254 x lo3 2 0.856 x lo2 0.645 x lo2 0.210 x lo2 0.834 x lo2 3 0.465 x lo2 0.319 x lo2 0.146 x lo2 0.457 x lo2 30 0.868 x 10 0.466 x 10 0.402 x 10 0.877 x 10 20 0.701 x 10 0.195 x 10 0.506 x 10 0.720 x 10 30 0.943 x 10 0.141 x 10 0.802 x 10 0.940 x 10 40 0.103 x lo2 0.128 x 10 0.899 x 10 0 . 1 1 0 ~ lo2 vector tensor k analytic diagonal L no constraint analytic diagonal L no constraint 1 0.213 x lop4 0.218 x lop4 0.235 x lop4 0.696 x lop4 0.699 x 0.788 x 2 0.179 x 0.182 x lop3 0.229 x 0.450 x 0.449 x 0.591 x lop3 3 0.107 x lop2 0.108 x 0.152~ 0.176 x lop2 0.178 x lop2 0.258 x 10 0.772 x lo-' 0.751 x lo-' 0.146 0.119 0.119 0.228 20 0.298 0.317 0.1 19 X 10' 0.490 0.497 0.173 x 10' 30 0.517 0.525 0.331 x 10' 0.912 0.896 0.592 x 10' 40 0.894 0.889 0.613 x 10' 0.151 x 10' 0.149 x 10' 0 .1 1 2 ~ lo2 the time correlation function C ( t ) = ( y k ( t ) Y k ( 0 ) ) e vs. t, where y k = q i q i and k = 5. Results from the various analytical approximations are shown along with results from the dynamical simulation. The value of C(0) must always be given correctly by any of the analytical approximations and is not displayed. The initial time derivative of the time correlations must always be given correctly by the full constraint formalism and the agreement seen in fig.1 for short times indicates only the (apparent) absence of programming errors. Of greater import is first that there is no statistically significant discrepancy between the constraint formalism and the dynamical simulations throughout the duration of the trajectories. The quantitative error caused by omission of constraints increases with mode number k, but it is significant even for k = 1 and for vector as well as tensor modes, unless L J p >> 1. Secondly, the mixing between modes demonstrated by the non-exponential decay is accounted for almost entirely by the equilibrium correlation To. Off-diagonal elements of L play essentially no role in the mixing of Rouse modes. According to much other evidence similar to fig.1, the agreement between analytical results and dynamical simulations seems to be satisfactory for constrained chains. We now proceed to examine various features of the analytical formalism and its consequences on the assumption that the constraint formalism is essentially exact for wormlike chains. Internal Friction It follows from the adequacy of the diagonal approximation to L, which will be further illustrated below, that the effect of constraints on the dynamics can be duplicated by an appropriate choice of a diagonal internal friction matrix A,Y in eqn (16). The formulae for Lkk (tensor modes) and L i k are given in eqn (17) and (18). The values of the required matrix elements A & are illustrated for tensor modes in table 1 and in fig.2. The analytical206 Brownian Dynamics of Chain Polymers 30 20 cz" 10 i ' " " " ' 0 0.2 0.4 0.6 0.8 k l N 3 Fig. 2. ( a ) The total mode friction constant Pk and ( b ) h X k , i.e. the contribution to Pk from constraints. Pre-averaged HI is included. Solid curves refer to 80 bond chains with L, == 48 and the dashed curves to 20 bond chains with L , / p =: 12. For each chain b = 5 , d s = 3, the angles are constrained at 122", and p = 8. results for the mode friction constants Pk = (BZk)-' are constructed from equilibrium simulations of L and To according to eqn (16)-( 18): P k B:; + h s k 2 k T r g k k / Lkj,. (19) h & P k - B;k1 2 k T r g k k / L k k - BZF' (20) The contribution to P k from constraints defines the internal friction h&: where I'; is the proper equilbrium correlation matrix of vector modes for the model, L is constructed for the required tensor modes with the appropriate constrained diffusion matrix B, and B,4 is constructed with omission of both constraints and internal friction. measures the effect of constraints in slowing down the relaxation; equivalently, A & is the amount of internal friction that has to be added to B:il, according to eqn (16), in order to make Bq reproduce the effects of constraints. Table 1 shows by the agreement between the exact P k and the approximate PL computed in the factored L approximation of eqn (14), that the factored approximation is reasonably good even for a very stiff chain.The approximation improves with increasing flexibility but it is reasonably good, as is the diagonal L approximation, even in the limit of a straight rod.P k is defined in such a way as to become proportional to the relaxation times T k for chains sufficiently long and flexible that To is diagonal as well as L. The relaxation times may be calculated generally, according to eqn (7), as the eigenvalues of ToL-', and a sampling of rates 1 / ~ ~ computed in this way is given in table 1. With the interpretation of A & just given, eqn (20) may be regarded as a calculation of the internal friction A & as a function of k, and the result can be contrasted with based on physical models or comparisons with experimental results. For long chains with elastic or constrained bond angles, h& becomes a function of k / NM.Fixman 207 as N + m . As the persistence length of the chain is increased, table 1 shows that h& for large k increases with k if bond angles are constrained (but it does not increase if the angles are elastic). For either model, elastic or constrained bond angles, h tjk increases for small k as the persistence length increases; see table 1. However, A & does not show the divergence as an inverse power of k / N , which was proposed in some of the works of Cerf, Peterlin and Allegra." Rather A & remains finite in the limit N -+ a. In this respect the results agree more closely with those of BazGa and Williams," who arrived at the simplest model in which htjk is independent of k and N. It seems probable that the rise in h& at small k has only an accidental relation to previous suggestions of this effect.Indeed, it also seems probable that a small change is required in the common picture of 'internal friction', or at least in the language used to describe it. In conventional language internal friction is introduced in order to slow down relative, internal motion in a real chain. It would follow that internal friction cannot affect the rigid-body translational or rotational motion of any chain. However, A & shows all the mathematical characteristics of internal friction and yet it is quite large for a rodlike molecule, which is capable only of rigid body motion. To dispose of this verbal paradox one must recall that internal friction is always used in practice as a correction to the motion of a 'Gaussian' model.This may be a generalized Gaussian model which makes the covariance matrix of bond vectors agree with that for a real chain or a rigid rod, but it is nevertheless a model which does not recognize dynamical constraints, i.e. the model does not accurately describe the velocities of the beads. This failure will apply also to the 'modes' of the Gaussian model, which are just linear combinations of bead coordinates. The velocity or relaxation rate of such a mode will not be calculated exactly in the Gaussian model even if the mode is a rigid body displacement or rotation. Discrepancies between velocities of the model chain and the real chain may emerge virtually instantaneously, because of constraints or vibrational forces, or after some time lag, because of hindered rotation barriers.In either case any initial displacement of the Gaussian chain will generate opposing forces different from those of the real chain and the relaxation rates of the Gaussian chain will generally be larger than for the real chain. For this reason 'internal friction' may be required to slow the relaxation (with the exception of a translational displacement without HI, for which the opposing force remains forever independent of configuration). First Cumulant and Translational Diffusion The first cumulant is defined" as R(q) = (p*Zp)J(p*p},, where N p = C exp (iq- &). k - 0 n(q) is an 'initial' time derivative of an observed structure factor, but the meaning of 'initial' depends on the degree of coarse-graining in time that is implied by the construc- tion of the diffusion operator 2, i.e.on whether or not constraints are used. We have performed the equilibrium simulations required for an evaluation of R(q) vs. q for various models with fluctuating HI. Results for constrained angle chains are similar to those found by Akcasu et al.," who considered freely jointed chains and used an analytical pre-average of M - ' . Results for a Gaussian model in which the bond vectors were generated independently, but the resultant bond lengths and angles were taken to be dynamically constrained, were litt!e different from those for a constrained 90" bond-angle chain. The unconstrained Gaussian model smooths out the peaks in the first cumulant relative to the unconstrained 90" chain, but has little effect on the location of the curve.Introduction of constraints into either the Gaussian or 90" chain smooths out the peaks, as does strong us. weak HI, depresses O ( q ) at large q, and makes the two models look quite similar.208 Brownian Dynamics of Chain Polymers Table 2. Translational diffusion constants as ratios of value with constraints to values without constraints. Constrained bond angles of 90” were used, constrained bond lengths b = 5 , and fluctuating HI. The Stokes diameters are shown as values of d s / b . The number of beads ranges from 14 to 112. The ratios have standard deviations of CCI. 0.003 N + l d S / b 14 28 56 112 ~ ~~ 0.5 0.977 0.966 0.965 0.961 1 .o 0.947 0.93 1 0.93 1 0.927 2.0 0.923 0.890 0.883 0.880 The depression in s l ( q ) for large q, owing to the addition of bond angle constraints, is rather modest for the relatively short chains that were studied.This is due to a significant contribution from the translational mode that would not be present for flexible chains of realistic size. Moreover, barriers to internal rotation would make the internal friction much larger than has been found for freely roating chains, and the high-frequency limit of is likely to be quite depressed. We estimate the depression for a simple model of internal friction. See also Allegra et al.19 We assume that N is sufficiently large that the sums over modes used to compute R(q) can be converted to integrals, that hydro- dynamic interaction is relatively unimportant for large Rouse mode numbers, and that A & = h for all k.This gives where Pd is here the friction constant of a single bead. If h = O the high-q limit of q-*sl( q ) is the value for Gaussian chains without constraints. With increasing internal friction the limiting value slowly decreases. The simulation results were typically about half that given by eqn (22) because of a substantial contribution from the translational mode k = 0 for short chains. Additional data on R(q) were gathered for rather large bead friction constants in order to assess the effect on the translational diffusion constant of coupling between the constraints and fluctuating HI. The full analytical constraint formalism was used for this purpose. It was anticipated that the constraints would behave similarly to internal friction2’ in affecting the diffusion constant. Similar results were indeed found, as is shown in table 2.The translational diffusion constant was considered to be proportional to R(O), and this amounts to using the Kirkwood formula as before2’ ( i e . the diffusion constant is taken p_roportional to a double sum of the averaged elements of the diffusion matrix), with G replacing G and with internal friction omitted. The diffusion constant so calculated is actually an upper bound to the true value.20 It seems reasonable on the basis of table 2 to conclude again that a coupling between fluctuating HI and local motion affects the translational diffusion constant even in the limit of infinite chain length. Stress Relaxation In fig. 3 and 4 the storage and loss moduli G’(o) and G”(o) are displayed as functions of the reduced frequency oR = wqo[ 77]/kT, where qo is the solvent viscosity and [ 771 is the intrinsic viscosity defined with concentration units of molecule volume-’.These modulii are designed [G’], and [G”]. in the book by Ferry.21 The more interestingM. Fixrnan 209 Fig. 3. Comparison of analytic approximations for the shear modulii us. wR. ( a ) The loss modulus G' and (6) the storage modulus G". Pre-averaged HI is included. The 40 bond chain has constrained bond angles of 163", 6 = 5 , and d s = 3. L,/p = 2. (-) Full analytic results; (- - -) constraint matrices omitted. The diagonal L approximation coincides with the full analytical result. -1.0 0.0 1.0 2 .o 3.0 log O R Fig. 4. Dependence of shear modulii on wR and chain length.Same as fig. 3 except that bond angle is 122" and only full analytic results are shown. (-) 80 bond chain with L,/p = 48; (- - -) 10 bond chain with L , / p = 6.210 Brownian Dynamics of Chain Polymers extreme approximations are shown along with the full analytical results in fig. 3. The complete neglect of constraints is obviously very unsatisfactory, particularly for the loss modulus, because the high-frequency limit of the viscosity vanishes in this approximation. The diagonal L approximation is generally quite satisfactory, although in some instances deviations are seen in G' at high frequencies, around wR= lo3. The storage modulus G'( w ) may become quite sensitive in the high-frequency region, around w 2 (1) in the reduced units, to sampling errors and approximations such as a restriction of L to its diagonal elements.In this high-frequency region, corresponding to t O( l ) , G ' ( w ) is levelling off at its asymptotic value. This value is determined by the residual flexibility at high frequencies, which is in turn very sensitive to coupling between Rouse modes under some conditions. These conditions are the ones which cause the constraints to leave very little freedom available to the local motions, as happens if the bond angles are constrained rather than elastic. Stress relaxation is the most sensitive of quantities ordinarily studied at high frequencies owing to its direct dependence on forces as well as configuration. The neglect of memory integrals in the constraint formalism is safest at the shortest times, but the numerical approximations and the physical applicability of the formalism become doubtful in an application to very short times.Ordinarily these short times will not be observed in mechanical relaxation experiments. Concluding Re marks The internal friction constant A & for mode k becomes a function of ( k / N ) rather than k and N separately, for moderate chain lengths with L,> p . This means that studies of sort chains suffice for the prediction of the properties of long flexible chains. For constrained but freely rotating chains the internal friction constants may be accurately evaluated from equilibrium simulations. In general the functional dependence of h& on k / N for constrained chains is different from previous suggestions.The function is not singular at small k / N, contrary to the behaviour of the total mode friction P k , which diverges as k-3'2 for chains with hydrodynamic interaction, and contrary to some of the suggestions of k-' or k-"2 divergence in A&. On the other hand, A & may show a sharp rise at small k / N , which is significant in the quantitative interpretation of stress relaxation in finite length chains, and a smooth rise at large k / N , which is significant to other experiments. Thus the most elementary physical models of internal friction seem inadequate. There remains the question of the effects of rotational barriers on these simple findings. The previous ~ i m u l a t i o n s ~ ~ of 90" constrained angle chains led to the pre- liminary suggestion that the Rouse-Zimm mode structure was well preseved in the presence of barriers.If this conjecture stands up to more thorough tests on stiffer chains than were studied before, as we expect, then the imitation of realistic backbone potentials with internal friction will be fully justified. However, one modification of the present treatment must be anticipated. The internal friction force will have to be written as an integral over the history of relative velocities, and the kernel will require a finite rather than an infinitesimal relaxation time. This work was supported in part by NIH GM27945, and was initiated as a result of conversations with R. Pecora and D. Roitman at a conference on the Rotational Dynamics of Small and Macromolecules in Liquids, at the University of Bielefeld, April, 1986, organized by Professors Th. Dorfmiiller and R. Pecora. References 1 P. E. Rouse Jr, J. Chem. Phys., 1953, 21, 1272; B. H. Zimm, J. Chem. Phys., 1955, 24, 269. 2 H. Yamakawa, Modern Theory of Polymer Solurions (Harper and Row, New York, 1971).M. Fixman 21 1 3 R. Zwanzig, J. Chem. Phys., 1974,60, 2717; M. Bixon and R. Zwanzig, J. Chem. Phys., 1978, 68, 1896. 4 M. Fixman and J. Kovac, J. Chem. Phys., 1974, 61, 4939. 5 D. B. Roitman and B. H. Zimm, J. Chem. Phys., 1984, 81, 6333; 6348. 6 K. Nagasaka and H. Yamakawa, J. Chem. Phys., 1985,83, 6480. 7 G. T. Evans, J. Chem. Phys., 1979, 70, 2362. 8 M. Fixman and J. Kovac, J. Chem. Phys., 1974,61,4950; M. Fixman and G. T. Evans, J. Chem. Phys., 1976, 64, 3474. 9 H. Yamakawa and T. Yoshizaki, J. Chem. Phys., 1981, 75, 1016; T. Yoshizaki and H. Yamakawa, J. Chem. Phys., 1986,84,4684; H. Yamakawa, T. Yoshizaki, and M. Fujii, J. Chem. Phys., 1986,84, 4693. 10 A. Z. Akcasu, B. Hammouda, W. H. Stockmayer and G. Tanaka, J. Chem. Phys., 1986, 85, 4734. 11 G. Allegra, J. Chem. Phys., 1986, 84, 5881. 12 Statistical Mechanics. Part B : Time -Dependent Processes, ed. B. Berne (Plenum Press, New York, 1977). 13 M. Fixman, unpublished results. 14 M. Fixman, J. Chem. Phys., 1978, 69, 1527; 1538. 15 M. Fixman, Macromolecules, 1986, 19, 1195; 1204. 16 M. Fixman, J. Chem. Phys., 1982, 76, 6124. 17 U. M. Titulaer and J. M. Deutch, J. Chem. Phys., 1975, 63, 4505. 18 E. R. Bazfia and M. C. Williams, J. Chem. Phys., 1973, 59, 2858. 19 G. Allegra, J. S. Higgins, F. Ganazzoli, E. Lucchelli and S. Bruckner, Macromolecules, 1984, 17, 1253. 20 M. Fixman, J. Chem. Phys., 1986, 84, 4085. 21 John D. Ferry, Viscoelastic Properties of’ Polymers (John Wiley, New York, 1980). Received 5th December, 1986