Stochastic Dynamics in Irreversible Nonequilibrium Environments. 2. A

thermosetting polymers.21-24 These systems play an important role in reaction injection ... Because of the self-similarity of polymer growth, macrosco...
1 downloads 0 Views 118KB Size
Download PDF
1070

J. Phys. Chem. B 1999, 103, 1070-1077

Stochastic Dynamics in Irreversible Nonequilibrium Environments. 2. A Model for Thermosetting Polymerization Rigoberto Hernandez* and Frank L. Somer, Jr. School of Chemistry and Biochemistry, Georgia Institute of Technology, Atlanta, Georgia 30332-0400 ReceiVed: September 3, 1998; In Final Form: December 2, 1998

A generalization of the generalized Langevin equation (GLE), the so-called irreversible GLE (iGLE) [Hernandez, R.; Somer, F. L. J. Phys. Chem. 1999, 103, 1064], is further extended to describe non-stationary environments in which the non-stationarity is induced by the macroscopic behavior of the ensemble itself, rather than an external force. Such a formalism lends itself to the dynamical study of the length distributions of growing polymers.

I. Introduction Flory1

Stockmayer2

The early statistical treatments of and provided a generally applicable description of the sequence distribution of polymers.3-7 Because living polymers,8-10 that is, polymerizing solutions in which the active sites remain unterminated, allow for equilibration of the sequence distribution, they have been a popular subject for theoretical investigations. Tobolsky and Eisenberg11 first treated equilibrium polymerization using mechanistic master equations. De Gennes12 and des Cloiseaux13 used renormalization group theory in interpreting continuum models of polymerization as a phase transition between small and high polymers. This interpretation was further validated by Wheeler and Pfeuty14 in showing that Scott’s generalization15 of the Tobolsky and Eisenberg model is equivalent to an Ising spin magnet in the limit that the spin vector dimension goes to 0. Several groups16-20 have studied equilibrium distributions and phase diagrams of living polymers by exploiting the isomorphism between continuum and lattice models. Another class of polymerizations in which one is interested in the dynamics of the sequence distribution is that of thermosetting polymers.21-24 These systems play an important role in reaction injection molding25,26 and have been the subject of large-scale finite-element calculations using semiempirical kinetic and visco-elastic equations (see, for example, ref 27). Because of the self-similarity of polymer growth, macroscopic kinetic equations provide a reasonably accurate reaction mechanism.28,29 What these studies often leave out is an understanding of the reaction dynamics when viscous, i.e., solvent, effects affect the reaction process. This is an operative mechanism in thermoset polymers where the polymerization ends not because the reactants have been depleted but because of diffusion effects due to the dramatic change in the viscosity of the fluid as a function of the extent of reaction. To achieve a better understanding of the thermosetting dynamics in which multiple branching leads to a highly crosslinked gel, in the present work we will consider the problem of chain polymerization in increasingly viscous environments. O’Shaughnessy and co-workers30,31 have constructed a FokkerPlanck master equation to describe the time evolution of the * Author to whom correspondence should be addressed. E-mail: [email protected].

sequence distribution. Their results exhibit the autoacceleration of polymer lengths, e.g., the Trommsdorff effect,32 which is characteristic of free-radical polymerization. The natural complement to this master equation is a stochastic model describing the dynamics of each member of the ensemble of growing polymers. The latter provides the rates and distributions, as a function of time as does the former, while also providing more detailed statistics such as the time-dependent autocorrelation function of individual polymer chain lengths in the presence of a polymerizing fluid. Unfortunately, the exact connection between the Fokker-Planck master equation and the Langevin equation for a nonstationarysnecessarily non-Markovians system is not currently known. In this work, we will show that the irreversible generalized Langevin equation (iGLE) presented in paper I33 of this series has the form that one would expect for such a Langevin equation, if the nonstationarity in the friction kernel is taken to be a function of the overall chain lengthening of the polymer ensemble. In future work, we will explore the precise connection between the proposed stochastic iGLE and the Fokker-Planck equation for polymerization. In the iGLE model for thermosetting polymerization, the friction of the bath is not merely space dependent34-38 and colored39,40 as in the usual GLE models.41 In addition, it depends on the average extent of reaction with some characteristic scaling exponent. Though this friction is implicitly colored, it cannot simply be subsumed under the time dependence of the stationary friction kernel because it will also depend on the conditions of the material at absolute times. As in the fruitful analogy seen in living polymers, thermosetting polymerization should be viewed as a phase transition which is driven by the reactivity of the polymer, and this is mimicked by the iGLE polymer model, in which the macroscopic growth self-consistently affects the detailed reaction dynamics. The potential of mean force (PMF) is constructed in section II in terms of an effective reaction path coordinate which is the contour length of the growing polymers. Time-dependent (nonstationary) solvent effects due to the change in the macroscopic properties of the polymerizing solvent are introduced in section III through the use of the iGLE. In section IV, the dynamics of this polymerizing model is analyzed in the chain-polymerization limit. Numerical results indicating the

10.1021/jp9836269 CCC: $18.00 © 1999 American Chemical Society Published on Web 02/03/1999

Irreversible Nonequilibrium Environments. 2

J. Phys. Chem. B, Vol. 103, No. 7, 1999 1071

structure of the iGLE with the nonstationary friction kernel, using biased potentials and the polymer PMF, are presented in section V. II. Potential of Mean Force In studies of polymers, the major emphasis has been on the understanding of the structure of the various materials given a particular size of the polymer, as determined by the number n of monomers in the polymer or n-mer. The focus has thus been placed on the free energy or potential of mean force for which the overall size R of an n-mer is unknown though its monomer length n is known, i.e.,

∫Ω dr e-βV(r)δ{R - |br n - br 1|}

e-βFn(R) ≡ Q-1

(1)

n

where β is 1/(kBT) for a temperature T, and the integral is over the space Ωn of all n-mers, without explicitly prohibiting monomer overlap, i.e., the phantom chain model. V is the interaction potential of the n-mers represented by the 3nr2, ... b rn), where b ri denotes the dimensional vector, r ≡ (r b1, b position of the ith monomer. Q is the partition function obtained by the integral without the δ-function constraint. For many models in which the bonds between the monomers are taken to be springs and in which long-range interactions are ignored, this free energy can be obtained exactly.3,4 The latter approximation entails ignoring the interaction between different polymers and even between distant monomers on the same polymer; in the spirit of the ideal gas law, the effect due to this interaction has come to be called the excluded-Volume effect.4,42 Scaling theories5,6,12 and renormalization-group techniques7,13,43-46 have added excluded-volume effects to Fn(R) and related averages, in obtaining deceptively simple scaling relations between polymer lengths and characteristic properties. A. Formal Expression. In the present work, we wish to shift the emphasis away from that of fixed n in order to allow for chemical reactions which extend n. The potential of mean force to be emphasized in this spirit is n-1

|b r i+1 - b r i|} ∑∫Ω dr e-βV(r)δ{R - ∑ i)1

e-βF′(R) ≡ Q′-1 ≡

(2)

n

n

∑n e-βF′(R) n

(3)

where R corresponds not to the size of the polymer but roughly to its contour length, and Q′ is the partition function of the monomer. The choice of Q′ sets the zero of free energy to be at R near l where l is the average monomer length. Equation 3 also serves to define the partial free energy F′n(R) associated with each n-mer space Ωn. The sum over n in the potential of mean force F′(R) is not divergent because, for n significantly larger or smaller than R/l the interaction potential will be large (i.e., F′n(R) will be large for R far from nl) and will thus provide exponentially small contributions to the integral. The PMF is smooth in R for nonrigid polymers because of the overall vibrations of the monomers, as well as the addition/ subtraction of a monomer. As long as the vibrations between the monomers are sufficiently tight that for a given R, the sum over n is dominated by n near R/l, we obtain the approximate equality,

F′(nl + ) ≈ F′n(nl + )

(4)

where || , l. This suggests that for R near nl the PMF is in a

Figure 1. The shape of the potential of mean force, F′n(R) is displayed as the solid curve. The biased potential (dashed curve) and the biased soft-core potential (dotted curve) used in section V are also displayed.

bound well supporting vibrational oscillations. For R between wells, the PMF will be dominated by the contributions from the two nearest polymer lengths,

e-βF′((n+σ )l+) ≈ (e-βF′n((n+σ )l+) + e-βF′n+1((n+σ )l+)) (5) †

†

†

where 0 < σ† < 1, and σ† is chosen so that (n + σ†) represents the transition states. Thus, one expects that the form of F′(R) is that of a series of wells centered near R ) nl, for integer n, with barriers in between. The self-similarity in the PMF, F′(R), is also necessary for this formalism to be in agreement with Flory’s assumption of equal reactivity.3 B. Phenomenological Expression. A direct calculation of F′(R) will be reserved for future work. The focus of the present paper is the use of a phenomenological model for F′(R) to explore its utility and applicability in describing polymer growth. For simplicity, a minimal model will be constructed which depends only on the characteristic barrier heights in the forward and backward directions, E†+ and E†-, respectively. Each F′n(R) can be approximated by a harmonic oscillator with a temperature-dependent frequency which is independent of n,

1 F′n(R) ≈ -n∆E + ω2(β)(R - nl)2 2

(6)

where ∆E (≡ E†--E†+) is the difference in energy between two neighboring wells. The sum in eq 3 will thus lead to a corrugated structure in F′n(R) with a series of wells at nl and separated by temperature-dependent barriers. As the temperature goes to infinity (β f 0), the corrugated structure will be lost both because ω(β) will go to zero and because the partial free energies will all contribute equally. At lower temperatures, however, the structure will not be lost, and for simplicity at a given temperature β-1, we will take the structure of F′n(R) to be that of a series of merged harmonic wells and barriers biased by ∆E, and satisfying continuity and differentiability. The precise form of this potential is presented in section VC and is displayed in Figure 1. It is parameterized by the forward and backward barrier heights, which are the only pieces of information needed to specify the potential interactions in the phenomenological polymer PMF model. III. Solvent Effects: iGLE The PMF with respect to an effective polymerizing (reactionpath) coordinate R′ as developed in the previous section provides

1072 J. Phys. Chem. B, Vol. 103, No. 7, 1999

Hernandez and Somer, Jr.

the average force on R′, but leaves out the solvent response of the nearby monomeric units. These neighboring monomeric units may be members of other chains (intermolecular), members of the same chain (intramolecular), or monomers in the solution (intermolecular). The need to include the intramolecular component is a consequence of the fact that we have assumed phantom chains. In principle, such solvent effects could be included with a generalized Langevin equation (GLE) which includes both the random kicks from the solvent as well as the frictional response. However, such a model does not account for the fact that, as the polymerization progresses, there is a macroscopic change in the composition of the material which leads to an increase in the friction. The irreversible GLE (iGLE) introduced in paper I33 provides us with a representation in which the friction can change according to a prescribed function g(t). The stochastic equation of motion in the iGLE in terms of R (with the mass set to 1, for simplicity) may be written as

R¨ +

dF′(R) + dR

∫t dt′g(t)g(t′)γ0(t-t′)R˙ (t′) ) g(t)ξ0(t)

(7)

where γ0 is a stationary friction kernel related to the random force ξ0 through the fluctuation-dissipation relation:

〈ξ0(t)ξ0(t′)〉 ) kBTγ0(t-t′)

(8)

As in the examples in paper 1, g(t) can be chosen as a switching function that increases the friction from a low equilibrium value characteristic of the unreacted material to a high equilibrium value characteristic of the final product. Such a function also obviates the need for a termination step in the model since the higher friction will lead to quenching in the polymerization. This is consistent with experimental observations in which the thermosetting process ends at curing (beyond gelation), well below 100% conversion. This friction is not the same as that which is usually discussed in the continuum chain models47-49 because here the friction corresponds to that along the effective reaction coordinate, which includes both small chain vibrations and monomer addition, and not between the polymers. Nonetheless, both of these are affected by solvent effects, and one presumes that there should be some connection between them. There is, however, a problem with choosing g(t) as a switching function: before the ensemble dynamics are run we do not know how the material has grown and hence what the profile of its characteristic macroscopic friction, vis-a-vis g(t), should be. Recalling that the viscosity of the material is related to the characteristic size of the polymers through a scaling relation,3 we posit that g(t) should be related to the characteristic length of the polymers through a similar scaling relation, i.e.,

g(t) ) |〈R(t)〉|ζ

(9)

where ζ is a scaling exponent. This choice of g(t) completely specifies the dynamics of the iGLE in eq 7. [Such a choice for g(t) does place a severe restriction on the behavior of the corresponding microscopic model to which the iGLE is a projection. However, as in similar remarks in paper 1, the use of g(t) captures much of the dynamical behavior that one expects from this system, without the additional complications imposed by more general models which are reserved for future work.] g(t) will behave like a switching function as long as 〈R(t)〉 quenches at long time. The latter must be true because eventually the growth of 〈R(t)〉 will lead to a large enough friction that the solvent response will quench; this will be illustrated in sections IV and V.

Notice that the choice of g(t) in eq 9 extends the analysis in paper 1 because in that discussion g(t) is assumed to be due to an external irreversible change in the solvent. Here, the collective motion of the ensemble effects g(t). Consequently, g(t) is self-consistently related to the dynamics. This may lead to concern that the arguments of paper 1 describing the relation between the nonstationary friction kernel and the random response may no longer apply. The following heuristic proof should allay such concerns: Once g(t) is known, by way of knowing 〈R(t)〉, then one can construct an iGLE which uses this g(t) as the external irreversible change in the solvent. By the arguments of paper 1, this latter iGLE must satisfy the connection between the nonstationary friction kernel and the random response. But since there is no difference between this latter iGLE and the iGLE with the self-consistent g(t), the iGLE defined with g(t) in eqs 7 and 9 also satisfies the appropriate equilibrium limits demanded of the iGLE constructed in paper 1. IV. Single-Hop Dynamics In one possible limit of the iGLE dynamics, the effective polymer ensemble will grow by successive single-addition reactions, i.e., chain polymerization. In GLE dynamics, this is referred to as single hops, in contrast with correlated hops50-52 in which the particles may hop over multiple barriers between quenches. In the present section the iGLE dynamics is restricted to the single-hop limit in order that we may provide a simple analytical (but approximate!) treatment to characterize the behavior. In section V, exact numerical calculations will be presented which illustrate the general behavior in which both single and correlated hops can both play a role. If we further assume that the dynamics are in quasiequilibrium as defined in paper 1,33 then at a given time t the average length of the ensemble may be written as

〈R(t)〉 ≡ n(th)l

(10)

where l is the average monomer length, ht is some characteristic time of the quasi-equilibrium interval Iht , and n(t) (≡ 〈R(t)〉/l) is the effective number of monomers at time t. Within this approximation, the friction kernels may be written locally as

γR(t,t′) ≈ n(th)2ζl2ζγ0(t-t′)

(11)

where the interval Iht contains both t and t′. In subsection IVA, the steady-state flux of the ensemble in quasi-equilibrium will be obtained for a given fixed n(th). If the system is always in quasi-equilibrium, and each of the associated intervals is long enough that steady-state conditions may be achieved, then the results in subsection IVA can be combined to obtain the flux behavior over all times ht . In the limit that the scaling exponent ζ is zero, the flux is independent of n(th), and so in subsection IVB, we obtain the propagation step characteristic of polymerization in time-independent environments. In general, the polymer-growth iGLE with arbitrary ζ will not exhibit steadystate conditions and quasi-equilibrium throughout the dynamics. Nonetheless, the early and late time behavior, in which the change in the friction is slow, does exhibit such conditions. In subsection IVC, the analysis of subsection IVA for these limits is used to characterize the overall dynamical behavior of the iGLE dynamics for polymer growth. A. Steady-State Flux. The mechanistic equation describing the rate process of the single-hop dynamics may be written for n > 0 as

Irreversible Nonequilibrium Environments. 2

J. Phys. Chem. B, Vol. 103, No. 7, 1999 1073

k+(th)

Mn {\ } Mn+1 k (t ) -

(12)

h

and the frequency λ† satisfies the transcendental equation

ω2b

†2

where Mn is a polymer with n monomer units M, and k+(th) [k-(th) is the rate of the forward [backward] reaction in the presence of the stochastic friction kernel during the time interval Iht . In the most general case of the dynamics, the interval will be too short to achieve steady-state like behavior, and the problem requires a treatment which includes the full timedependence in the reaction rates. In the present section, we treat the simple limit in which Iht is long enough to achieve steadystate. In this limit, the problem can then be solved by combining the solutions for each such time interval. The corresponding master equation within a given steady-state interval may be written as

P˙ 1 ) -k+P1 + k-P2

(13a)

P˙ n ) k+Pn-1 - (k+ + k-)Pn + k-Pn+1, for n > 1

(13b)

where Pn(t) is the population (concentration) of Mn, and the dependence on ht has been dropped for notational simplicity. Jung and Berne52 studied the multibarrier crossing problem with an external dc field. Their problem is identical to that in eq 13b, if we ignore the terms in their master equation which correspond to correlated hops, and we ignore the edge effect due to P1 in eq 13a. The kth eigensolution of eq 13 can be formally written as -λkt (k) cn P(k) n (t) ) e

(14)

where λk and c(k) are the corresponding eigenvalue and eigenvector of the linear matrix contained in eq 13, respectively. Either through a direct calculation of the flux using these eigenvectors, as was done by Jung and Berne,52 or by inspection of eq 12, the steady-state flux at all the barriers is

j ) 2π(k+ - k-)

(15)

which is independent of n > 1. The edge effect at n ) 1 can be resolved by requiring the boundary condition that the flux at n ) 1 be equal to the steady-state flux in eq 15. This is tantamount to the usual requirement in activated dynamics: reactants are fed into the initial well to maintain the initial population. In order to calculate the flux, we need to obtain the rates, k+ and k- in the presence of the friction kernel of eq 11. Within the quasi-equilibrium regime in which eq 10 holds, these rates are precisely the rates of the Kramers turnover problem,53 which may be written as39,54,55

( ) ( )

ω0 -βE†+ e k+(th) ) κ+(th) 2π

(16a)

ω0 -βE†e k-(th) ) κ-(th) 2π

(16b)

where E†+ (E†-) is the activation barrier for the forward (backward) reaction, and ω0 is the harmonic frequency of the well. The correction to transition state theory, κ(, contains the friction dependence, and may be written as56,57

[∫

1 λ† κ((th) ) exp ωb π

∞

]

2 dy ln{1 - eδ((1+y )/4} 2 1+y

-∞

(17)

where iωb is the imaginary harmonic frequency of the barrier,

λ )

(18)

1 + n(th)2ζl2ζγ0[λ†]/λ†

The friction term also enters in kBTδ( which is the average energy loss of the subsystem between visits to the barrier.57 It will not be the same in the forward and backward direction, as the forward and backward barrier heights are not equal. The flux for these systems can thus be written as

j ) ω0[κ+(th)e-βE+ - κ-(th)e-βE- ] †

†

(19)

The terms κ( turn over from small values at low friction, to approximately unity at intermediate friction, and finally to small values at high friction. Since the friction is proportional to n(th), the turnover in the rates will proceed according to n(th). Besides providing a justification of the claim that the proposed model can mimic the dynamical behavior of the polymerization process, the present analysis also provides estimates of the rates which will be used in future work to characterize the dynamics. B. Simple Kinetics. In the limit that ζ ) 0 there is no dependence on the polymer size, and the fluxes in eq 19 are independent of time. Moreover, if E†- . E†+, the forward reaction will dominate. This is tantamount to rewriting the mechanistic equation as the propagation step k0

Mn 98 Mn+1

(20)

where the flux in eq 19 can now be taken with only the first term in the RHS, i.e., the forward reaction only. The flux is always constant and positive, leading to infinite propagation (polymerization lengths) unless a termination mechanism is added to the master equation. Thus we have recovered the usual propagation step of polymerization reactions. C. Autoacceleration. In the present work, however, a termination step is not needed because the increasing friction quenches (terminates) the reaction process. This is the regime in which the scaling exponent leads to an increase in the friction with n, i.e., ζ > 0. If the early-time friction kernel γR(0,0) is also taken to be small, then at early times with n near 1, the dynamics will be energy-diffusion limited and consequently slow. Within this pseudo-steady-state, the flux expression of eq 19 applies near the peak of the length distribution. The flux is positive and will lead to increased n with time. This, in turn, leads to an increase in the flux, exhibiting the aforementioned Trommsdorff effect32 in which polymerization autocatalyzes its own growth. Eventually, the rates will be large enough that the steady-state picture will no longer apply. Nonetheless the dynamics will continue to grow a larger polymer ensemble until the fluxes once again depress to the pseudo-steady-state regime in which eq 19 once again applies near the mean of the distribution. This is now the diffusion-limited regime, and the fluxes will continue to depress with increasing polymerization. This suggests that the polymerization in this picture never fully ends, but rather the polymerization simply slows down to such slow rates that no appreciable growth is observed, i.e., it quenches to a glassy state. While this analysis explicitly applies to chain polymerization (because of the requirement of single hops), the dynamics of step polymerization (including correlated hops) will be exhibited in the following section, with the qualitative behavior remaining the same. V. Numerical Results and Discussion In this section, we will show that the iGLE, with the friction kernel via g(t) in eq 9 and the polymer PMF, does exhibit

1074 J. Phys. Chem. B, Vol. 103, No. 7, 1999

Hernandez and Somer, Jr.

dynamics characteristic of chain and step polymerization as described in the previous section. The former are represented within the dynamics as single hops whereas the latter are represented as correlated hops. The numerical details regarding the integration of the iGLE with an external g(t) can be found in paper I;33 the additional details necessary to integrate the iGLE with a self-consistent [vis-a-vis g(t) ] friction kernel are described in section VA. The validity of this nonstationary friction kernel is illustrated in section VB with an unstructured biased potential (constant force). Using a simple model for the polymer PMF, the iGLE polymer dynamics is illustrated in section VC. A. Integrating the iGLE. The presence of the 〈R(t)〉dependent friction in the iGLE complicates the numerical procedure because 〈R(t)〉 is not known a priori, and thus the stochastic dynamics must be solved self-consistently. Because only 〈R(t)〉 for t less than the present times are required, however, one can adopt one of the following numerical procedures: 1. A fully self-consistent iteration scheme: (a) use a trial function for 〈R(t)〉, (b) run the stochastic dynamics with the trial function for several particles, computing a new 〈R(t)〉, (c) check if the new 〈R(t)〉 agrees with the old, (d) if not, use the new 〈R(t)〉 as the trial function in step (b), and iterate until convergence. 2. An on-the-fly scheme in which 〈R(t)〉 is calculated for a large swarm of particles during the dynamics and 〈R(t)〉 for the previous known times is used to compute the future dynamics for each effective particle. 3. A hybrid scheme in which the current 〈R(t)〉 for the swarm of particles in the self-consistent iteration scheme is averaged with the trial function in order to speed up the convergence. The self-consistent procedure has the advantage that it will have almost perfect scaling with the number of particles and is thus extremely parallelizable. The on-the-fly scheme requires only one time-dependent nonparallel run, but it is limited in accuracy by the number of trajectories that can be run simultaneously, which is in turn limited by the memory size of the computer. The hybrid approach affords a compromise between memory and computing resources. In the present calculations, only the on-the-fly procedure was used and the runs took on the order of hours to days on Silicon Graphics R10000 processors. B. iGLE with Biased Potential. In order to explore the dynamics of the iGLE with the 〈R(t)〉-growth-dependent nonstationary friction vis-a-vis g(t) in eq 9, we first take the potential to be that of the biased potential

V(R) ) -fbR

(21)

where fb is a constant applied force, and the domain of R is over all space. The strength of the applied force fb is taken to be 5kBT/l where l is some characteristic size for the system. As in the polymer growth model to be presented in section VC, l will be taken to be characteristic of a monomer size, and R will be written in dimensionless units of l. As in paper 1,33 the simulation parameters are chosen to be N ) 10 000, kBT ) 2.0, γ0(0) ) 10.0, τ ) 0.5, ∆t ) 4.0 × 10-3, and ∆tξ ) 4.0 × 10-4, all in arbitrary units. (Recall that ∆t is the time step for the iGLE dynamics, and ∆tξ is the time step for the auxiliary dynamics of the associated Langevin equation with respect to γ0(t) ) γ0(0)exp(-t/τ).) The scaling exponent ζ in the friction kernel of eq 9 is set to 2. This value has been chosen to be small enough that the lengths seen in the polymer PMF

Figure 2. The average position 〈R(t)〉 is shown as a function of time for an on-the-fly iGLE dynamics with a swarm of N ) 10 000 particles localized at R ) 0 at t ) 0.

Figure 3. The friction kernels normalized to γ(t,t′) for the biased potential at various later times t are displayed as a function of the time difference t - t′. The displayed times t are 0.8 (dotted), 1.2 (dashed), 1.6 (long-dashed), 2.4 (dot-dashed), and 8.0 (solid). The corresponding values of γ(t,t) are 127.0, 36.5, 58.0, 102.1, and 378.2, respectively.

dynamics (section VC) are significantly longer than dimers and trimers while being large enough that the simulations converge quickly. The average length 〈R(t)〉 for the biased potential with the initial condition that all the particles are localized at the origin at time 0, is shown in Figure 2. Because the applied force fb is relatively large and the initial friction is 0, the initial behavior is essentially ballistic. It then oscillates to a damped diffusive growth. The normalized friction kernel at various times t of the simulation with the biased potential is displayed in Figure 3. These show dramatic deviations from the structure of the normalized stationary friction kernel γ0 at early times, and go smoothly toward that of the normalized γ0 at later times well before the distribution quenches. Note that, at t ) 8.0, the friction has reached the quasi-equilibrium limit in which γ(t,t′) ∝ γ0(t - t′) in the (t - t′) region that γ(t,t′) is far from zero, even though the static friction is still increasing with absolute time t. (The latter observation can be seen from the fact that 〈R(t)〉 in Figure 2 is increasing at long times and consequently the static friction in eq 9 is also increasing at long times.) In this regime, the response of the solvent due to the change in the static friction is changing slowly compared to the response time of the solvent in quasi-equilibrium.

Irreversible Nonequilibrium Environments. 2

J. Phys. Chem. B, Vol. 103, No. 7, 1999 1075

{

1 2 for R < Rm ω0(R - l)2 2 V′sc(R) ≡ -fb(R - l) - f2b/(2ω20) otherwise

(22)

where Rm (≡ l - fb/ω20) is the value at which the soft-core parabola and the line are matched, satisfying continuity and differentiability. ω0 is the well frequency and fb is the biased force as in the previous subsection. In the phenomenological polymer PMF used in the present models, we will take the barrier frequency ωb equal to the well frequency ω0. Noting the self-similar structure of the polymer PMF, the simplest model that still retains the barrier structure is that of a series of merged parabolas which is the natural extension of the double-well model of Straub et al.58 For completeness, it may be written as Figure 4. The velocity autocorrelation function for the biased potential is displayed as a function of time at t equal 4.0, 6.0, and 8.0.

{

V′PMF(R) ≡ 1 2 ω (R - nl)2 - fb(n - 1)l 2 0 1 for R < R′m + l with n ) 1, 2 1 or R′m + n - l < R e R′m + nl 2 1 - ω20(R - nl - 2R′m)2 - fb(n - 1)l + E†+, 2 1 for R′m + nl < R e R′m + n + l 2

(

)

(

(23)

)

where n (g 1) is implicitly obtained, and the first case applies to wells, and the second to barriers. The relative position of the transition state and the frequency are specified by

fb l R′m ) - 2 4 ω

(24a)

0

Figure 5. The mean-square velocity for the biased potential is displayed as a function of time, showing that it obeys the equipartition relation for intermediate to long times.

ω0 ) 2fb 2E†+ + fbE†+ - 2E†+ x1 + fbl -1/2

The change in the dynamical correlation is illustrated in Figure 4 where the velocity autocorrelation function is displayed at three different times. The correlation times are clearly changing with the increasing friction. Because of the strong deviations from stationarity, the equipartition value (2.0) of the mean-square velocity is not obtained throughout the dynamics as shown in Figure 5. At early times, before the friction kernel has increased to a value that is large enough to dampen the dynamics, one sees strong deviations from equipartition. However, at intermediate to long times, the friction grows strong enough to dampen the motion, and equipartition is established within this quasi-equilibrium regime. C. iGLE with Polymer Growth PMF. The dynamics of the polymer PMF within the iGLE can be illustrated with a simple phenomenological potential that is characterized only by the barrier height between the n-mer and the (n + 1)-mer and the exoergicity as characterized by the biasing force fb. Without these barriers the polymer PMF reduces to the biased potential of the previous subsection except for the fact that the relevant domain is now that of nonnegative R corresponding to finite polymer lengths. As such, a biased soft-core potential that is linear everywhere except near the core provides dynamics which do not enter into the nonnegative region, and serves as a better comparison to the dynamics with the polymer PMF. The soft core potential is taken to be

respectively. As a consequence, the only parameters of the phenomenological polymer PMF are the monomer size l (which is taken to be 1 for simplicity), fb (≡∆E/l), and E†+. The specific parameters for the iGLE simulations are N ) 100, kBT ) 2.0, fb ) 1, γ0(0) ) 1.0, τ ) 10, ∆t ) .007, ∆tξ ) .0007, and ζ ) 2. Note that fb is smaller here than in section VB, thereby deemphasizing the exoergicity. In the present simulations the barrier height E†+ in the forward direction is taken to be either 4kBT, 6kBT, or 8kBT. The normalized friction kernel experienced by the polymer PMF is displayed in Figure 6 at an early time. Because the decay time τ of the stationary friction γ0 is small compared to the timescales of the dynamics in 〈R(t)〉, the normalized friction kernels exhibit little structure relative to the normalized γ0(t-t′). However, the absolute friction constant is increasing during these runs, and this corresponds to a dynamics which is in quasi-equilibrium at short times (t