3456
J. Phys. Chem. B 2000, 104, 3456-3462
Stochastic Dynamics in Irreversible Nonequilibrium Environments. 4. Self-Consistent Coupling in Heterogeneous Environments Frank L. Somer, Jr. Department of Chemistry, St. John’s UniVersity, 8000 Utopia Parkway, Jamaica, New York 11430
Rigoberto Hernandez* School of Chemistry and Biochemistry, Georgia Institute of Technology, Atlanta, Georgia 30332-0400 ReceiVed: August 11, 1999; In Final Form: January 25, 2000
A generalization of the generalized Langevin equation (nonequilibrium dynamics) has recently been introduced in order to model chemical reactions which take place in environments with irreversible density and temperature changes (Hernandez, R., Somer, F. L. J. Phys. Chem. 1999, 103, 1064, 1070). The so-called irreversible generalized Langevin equation (iGLE) has been shown to be relevant to reaction polymerization when the irreversibility in the equilibrium solvent response is self-consistently coupled to the macroscopic average extent of polymerization, itself. In the present work, the iGLE is extended to allow the reacting polymers to be coupled to a finite number w of neighboring polymers, thereby allowing for heterogeneity in the environment response.
I. Introduction An irreversible generalized Langevin equation (iGLE) has been developed to describe reaction polymerization (or other processes occurring within irreversibly driven environments) in which the dynamics is quenched by a nonstationary increase in dissipation vis-a-vis the friction kernel.1-3 This increase may be self-consistent in the sense that it may depend on the collective growth of all the polymers through the average of the polymer lengths with a given scaling exponent.2 However, the reacting polymers are presumably not in a uniform polymeric environment. Instead, each is in a fluctuating local environment whose irreversible dissipative response is determined primarily by interactions with the finite number of effective polymers, w, in this local environment. The objective of this work is thus to develop a new version of the iGLE, of degree w (WiGLE), which can describe the dynamics of a tagged particle within a self-consistent heterogeneous irreversible environment. The self-consistent iGLE theory and a heuristic argument extending it to the WiGLE are presented in section II. Stochasticdynamics simulations for a variety of potentials of mean force (including one representative of reaction polymerization) are presented in section III, illustrating the utility of the WiGLE formalism as well as providing empirical support for its validity. One advantage of WiGLE, which is clearly seen in the results of section III, is that it lends itself to a much larger number of effective polymers than the iGLE, given the same amount of computing resources. It should also be noted that this is the last paper in this series1-3 in which we have constructed a progression of phenomenologicalsbut internally consistentsschemes for incorporating nonstationary environments into stochastic dynamics. In a separate paper,4 one of us has shown that the constanttemperature iGLE of refs 1 and 2 may be obtained as the projection of a microscopic Hamiltonian system, albeit a * Author to whom correspondence should be addressed. E-mail:
[email protected].
nonconservative one. A similar microscopic justification of the variable-temperature iGLE of ref 3 and the WiGLE of the present paper will be the subject of future work. II. iGLE to WiGLE A. Review of iGLE. The iGLE is a stochastic differential equation describing the motion of a chosen coordinate, which in one dimension may be written as
V˘(t) ) -
∫t dt′ γ(t, t′) V(t′) + ξ(t) + F(t)
(2.1)
where F(t) (≡-∇qV(q(t))) is the external force, V ()q˘ ) is the velocity, and q is the mass-weighted position along the chosen coordinate. Note that the chosen coordinate may represent the actual position of a chosen Brownian particle, the distance between two interacting particles, the motion along some effective reaction path of collective motions, etc. In the particular cases explored in the numerical simulations in section III, where q is interpreted as an effective reaction path for polymerization, q will be denoted as R. The nonstationary friction kernel γ(t,t′) and random force ξ(t) are related through a nonstationary generalization of the fluctuation-dissipation relation (FDR)
〈ξ(t)‚ξ(t′)〉 ) kBTγ(t,t′)
(2.2)
where kB is Boltzmann’s constant and T is the temperature of the bath. The nonstationarity in the iGLE is characterized by a time-dependent function g(t) such that
γ(t,t′) ) g(t) g(t′) γ0(t-t′)
(2.3)
where γ0(t-t′) is a stationary friction kernel related to a a random force ξ0(t) through the FDR
〈ξ0(t)‚ξ0(t′)〉 ) kBTγ0(t-t′)
(2.4)
For a given g(t), the iGLE construction is thus completely specified by requiring ξ(t) ) g(t) ξ0(t). In the limit that g(t)
10.1021/jp9928762 CCC: $19.00 © 2000 American Chemical Society Published on Web 03/22/2000
Irreversible Nonequilibrium Environments
J. Phys. Chem. B, Vol. 104, No. 15, 2000 3457
becomes constant and equal to gj, this relation also guarantees that the system will be described by an ordinary GLE with the friction kernel given by gj2γ0(t-t′). In previous work1,2 we have shown that the iGLE, interpreted as a nonstationary (“irreversible”) GLE, satisfies the correct equilibrium behavior in quasi-equilibrium limits as well as more generally in illustrative models. In the case that the solvent response changes due to outside forces, the form of g has been taken as a switching function which changes the solvent from a lower to a higher effective friction constant, γ0. As will be discussed below, we have also explored the case in which the change in solvent response is driven by changes in the internal state of the system, under both constant and biased potentials. These studies have graphically demonstrated that the iGLE dynamics satisfies equipartition well beyond the equilibrium limit. B. Extension to WiGLE. In the iGLE, each realization of the chosen coordinate’s motion is manifested within an equivalent environment, and consequently there is no need to differentiate these realizations. In the WiGLE, however, each of N realizations of the chosen coordinate labeled by n is in its own unique environment, and consequently must be determined by a unique stochastic differential equation
V˘n(t) ) -
∫t dt′ gn(t) gn(t′) γ0(t-t′) Vn(t′) +
gn(t) ξ0(t) + F(t) (2.5)
where n ranges from 1 to N. By analogy with the equations in the previous subsection, the friction kernel is connected to the random forces through a fluctuation-dissipation relation
〈ξn(t)‚ξn(t′)〉 ≡ kBTgn(t) gn(t′) γ0(t-t′)
(2.6a)
) kBTγn(t,t′)
(2.6b)
The primary new addition in the present formalism is that the nth chosen coordinate is retarded by a heterogeneous, irreversibly driven friction characterized by
gn(t) ≡ 〈|q(t)|〉ζn 〈|q(t)|〉n ≡
1
∑
w + 1 i∈Sw,n
(2.7a) |qi(t)|
coordinate q as the mass-weighted length of a polymer R. In earlier work, we have shown that R is a reasonable choice to represent the reaction path of polymerization.2 To fully solve the WiGLE for a nontrivial degree w, one needs to identify the neighbor sets Sw,n which characterize the motion of the effective reaction coordinate of the nth realization of the polymer, Rn(t). In a three-dimensional dynamical system this is not a simple task because the neighbor sets are themselves dynamic, and any one-dimensional ordering (through some given collective index) of the polymers will invariably result in contacts between distant indices. In a quasi-nematic twodimensional system, however, the polymers can be indexed by a simple collective index, with quasi-static neighbor lists Sw,n containing indices near n. Representing the corresponding coupling though a matrix, the latter would be described by a banded matrix with w off-diagonal stripes, whereas the former is an unstructured, but sparse, matrix. For illustrative purposes, most of the numerical simulations of this section involve a WiGLE with banded coupling. As we have described, this is one possible limit of the structure of the environment. It may also be appropriate for a general system in the sense that the dynamical effects are dominated by the degree w and not the specific makeup of the neighbor lists, Sw,n. To this end, we also present results for a blocked WiGLE in which the polymers are coupled only to those polymers in a fixed neighborhood of size w. A. WiGLE with Constant Potential. As in the previous papers in this series,1-3 we take the stationary part of the friction kernel to have the exponential form9-12
γ0(t-t′) ) γ0(0)e-|t-t′|/τ
and generate the random force Via the auxiliary Langevin equation
1 ξ˙ 0(t) ) - ξ0(t) + ξG(t) τ
III. Numerical Simulations In order to illustrate the rich dynamics of the WiGLE concretely, throughout this section, we specify the chosen
(3.2)
where ξG is a Gaussian random process characterized by the white-noise correlation function
(2.7b)
where ζ is a scaling exponent characteristic of the growth of the friction with the local environment dynamics, and Sw,n(3 n) is the set of labels of the w + 1 realizations in the local environment of the nth chosen coordinate (i.e., qn(t), itself, plus its w tagged neighbors). In the limit that w f N - 1, the average in eq 2.7b will go to the macroscopic average of q, and each of the the realizations of the chosen coordinate will be in the same effective irreversible environment. This limit is precisely what was studied earlier with self-consistent dynamics.2 In the limit that w ) 0, the WiGLE becomes an uncoupled stochastic differential equation which is self-dissipated in the sense that its own evolution will increase the friction of its environment. That is, when w ) 0, we recover the generalized Langevin equation with spacedependent friction.5-8 In between these limits, the iGLE represents a swarm of effective (projected) realizations, each of whose environmental response is affected by the overall motion of the neighboring realizations.
(3.1)
〈ξG(t) ξG(t′)〉 ) 2
( )
kBT γ (0) δ(t - t′) τ 0
(3.3)
A set of representative parameters have been used in all the reported calculations; specifically, these are N ) 513, kBT ) 2.0, γ0(0) ) 10.0, τ ) 0.5, ∆t ) 4.0 × 10-3, and ∆tξ ) 4.0 × 10-4, all in reduced units. [Recall that ∆t is the time step for the iGLE dynamics and ∆tξ is the time step for the integration of the auxiliary Langevin equation, eq 3.2.] The applied force F(t) is set to 0, and as in ref 2 the scaling exponent ζ in the friction kernel of eq 2.7a is set to 2. WiGLE simulations have been performed for various values of the correlation width ws 0, 4, 16, 64, and 512sand with initial conditions for the positions at Rn(0) ) 0 and the velocities Vn(0) chosen from a Maxwell-Boltzmann distribution. The mean square of the irreversible change in the environment, 〈gn2(t)〉, is displayed in Figure 1. The rapid growth at early times results from the nearly ballistic motion at very small values of 〈|R(t)|〉n. As the friction grows, the memory of the friction kernel gives rise to oscillations in 〈gn2(t)〉, the amplitude of which is attenuated with time. The amplitude and persistence of these oscillations decreases with w until at w ) 0 they are nearly obscured by thermal noise. This rapid loss of coherence is due to increased heterogeneity among the individual friction kernels.
3458 J. Phys. Chem. B, Vol. 104, No. 15, 2000
Somer and Hernandez
Figure 1. Mean square of the irreversible change in the environment, 〈γn(t,t)〉 ) γ0(0)〈gn2(t)〉, for the constant-potential case shown here as a function of time for various values of w.
Figure 3. Variance in the irreversible change in the environment, 〈γn2(t,t)〉 - 〈γn(t,t)〉2, for a constant potential displayed as a function of time at various values of w.
Figure 2. Average nonstationary friction kernel 〈γn(t,t′)〉 for the constant-potential case displayed at several different times t as a function of previous times t - t′, and at various values of w. In all cases, the abscissa is displayed as a ratio normalized to the corresponding value of 〈γn(t,t)〉 (which may be read off Figure 1), so as to easily display them in one figure, and the time t to which each line corresponds is indicated in the legend.
Figure 4. Velocity autocorrelation function, 〈V(t) V(t′)〉, for a constant potential displayed at each value of w and initial times t ) 4 (dashed curve), 8 (dotted curve), and 16 (solid curve), as a function of previous times t - t′. The curves are labeled by the corresponding value of w; and the t ) 8 and t ) 16 curves are offset by four and eight time units, respectively, for ease of viewing.
The average friction kernel, normalized to γ(t,t), is shown at various times in Figure 2. For all values of w we see considerable irreversibilitysas measured by the degree of divergence of γ(t,t) from the exponential form of γ0(t-t′)sat short times. The irreversible behavior appears to be more pronounced and longer lived the smaller the value of w. This is once again due to increased heterogeneity in the local particle environments: for small w, particles at small displacements R will encounter relatively low friction, thus allowing nearly ballistic motion that will enable large relatiVe changes in gn(t) over the time scale of the friction kernel, as is necessary to show substantial nonstationarity. This heterogeneity is further illustrated in Figure 3, where we show the time dependence of the variance of the local equal-time value of the friction kernel, γn(t,t). While the mean of γn(t,t) (which is proportional to g2(t)) appears to be fairly independent of w for the early times shown
in Figure 1, we see in Figure 3 that its variance increases dramatically with decreasing w. This is due to fact that gn (and therefore the corresponding friction kernel) differs from gn+1 only by the replacement of a single term in the sum in eq 2.7b: the greater the value of w, the less the relative importance of this single differing term, and hence the greater correlation between the friction kernels of neighboring particles. The velocity autocorrelation function (VAF) at three different times, for various values of w, is shown in Figure 4. At larger values of w, the frequency of the oscillations in the VAF increases with increasing friction, as a function of time. The oscillations fall off much faster as w decreases, until at w ) 4 they are swamped by noise after the initial minimum. The markedly slower response at w ) 0 is once again seen as a result of the motion of slower moving particles being unaffected by the positions of faster particles, such that they remain in
Irreversible Nonequilibrium Environments
J. Phys. Chem. B, Vol. 104, No. 15, 2000 3459
Figure 6. Mean square of the irreversible change in the environment, 〈γn(t,t)〉 ) γ0(0)〈gn2(t)〉, for the biased potential of eq 3.4 shown here as a function of time for various values of w.
Figure 5. Mean-square velocity, 〈V2(t)〉, for a constant potential displayed for various values of w as a function of time t.
low-friction environments longer. This situation also affects the evolution of the mean-square velocity (MSV), as shown in Figure 5. For the autocorrelated (w ) 0) case we see that at very short time the MSV drops rapidly from its equipartition value of 2, falling to nearly half that before returning to the equilibrium value around t ) 5. This initial falloff results from the rapid loss in kinetic energy of the initialized fast particles (which quickly move into high-friction regimes and consequently experience fast energy relaxation) and the slow increase in kinetic energy of the initialized slow particles (which remain in the low-friction environments for longer time and are consequently slow to reach thermalization). This effect is lessened as w becomes larger and makes each friction kernel less responsive to its corresponding particle’s trajectory, such that by the time we reach w ) 16 the initial dip in the MSV, while still present, is reduced to roughly the amplitude of thermal fluctuations. B. WiGLE with Biased Potential. In this section, we explore the behavior of the WiGLE under the influence of a constant applied force
F(t) ) fb ≡ 5kBT/l
(3.4)
where l is some characteristic size for the system. The other parameters and initial conditions are as in the previous section. The time dependence of 〈gn2(t)〉, shown in Figure 6, is qualitatively similar to that seen in the free-particle case. The growth is faster and the oscillations are more pronounced as a result of the applied force, but as in Figure 1 the amplitude of the oscillations decreases with w, because of the increase in the heterogeneity of the local environments. The biasing force also causes the average friction kernel, displayed in Figure 7, to show much more dramatic irreversible behavior at short timessthough, as before, less so for smaller values of wsthan was seen in the constant-potential case. However, since the particles are now forced into high-friction environments very quickly, regardless of initial velocity, the quasi-equilibrium regime is reached more quickly under the biased potential. The degree of heterogeneity in the local environments, as illustrated
Figure 7. Average nonstationary friction kernel 〈γn(t,t′)〉 for the biased potential is displayed at several different times t and various values of w as a function of the previous times t - t′. In all cases, the abscissa is displayed as a ratio normalized to the corresponding value of 〈γn(t,t)〉 (which may be read off Figure 6), so as to easily display them in one figure, and the time t to which each line corresponds is indicated in the legend.
in Figure 8, shows the same sort of dependence on w as seen for the free-particle case. In fact, apart from greater short-time oscillations for each value of w, the curves in Figure 8 are qualitatively similar to the corresponding curves in Figure 3. This suggests that the relative degree of heterogeneity might in fact be dominated by the magnitude of w, and not the applied force. The VAFs for the biased-potential case are shown in Figure 9. The anomalously slow response at w ) 0, seen for the freeparticle case in Figure 4, is absent here since all the particles are quickly forced into high-friction regimes. Otherwise, the behavior is qualitatively similarsthough more pronounceds with the correlation times falling dramatically with decreasing w. As seen in Figure 10, the biasing force, in conjunction with the low initial frictions, causes a rapid increase in the MSV at
3460 J. Phys. Chem. B, Vol. 104, No. 15, 2000
Somer and Hernandez
Figure 8. Variance in the irreversible change in the environment, 〈γn2(t,t)〉 - 〈γn(t,t)〉2, for the biased potential displayed as a function of time for various values of w.
Figure 10. Mean-square velocity, 〈V2(t)〉, for the biased-potential case displayed for various values of w as a function of time t.
Figure 9. Velocity autocorrelation function, 〈V(t) V(t′)〉, for the biased potential displayed at various initial times t and values of w as a function of the previous times t - t′. The curves are labeled by the corresponding value of w. The format of the plot is as in Figure 4, with the VAF at initial times t ) 4, 8, 16 displayed for each value of w.
early times. The magnitude of the initial peak decreases with w, again as a result of each friction kernel becoming more responsive to the motion of its corresponding particle. The subsequent memory-induced peaks in the MSV, clearly evident for the larger values of w, are smeared out for smaller values where the greater inhomogeneity among the local environments weakens the correlations between velocities of different particles. In all cases, equipartition is satisfied after a short equilibration time of t ∼ 5. C. WiGLE with Biased Soft-Core Potential. In this section we turn our attention to the biased soft-core potential introduced in ref 2:
{
1 2 for R < Rm w0 (R - l)2 V′sc ≡ 2 -fb(R -l) - fb2/(2ω02) otherwise
(3.5)
where Rm (≡l - fb/ω02) 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 section. As discussed there, this can be viewed as a simplified potential of mean force (PMF) for polymer growth in which the “polymers” grow without barriers. In addition to the banded WiGLE studied in the previous sections, here we also explore the block WiGLE in which each particle is influenced only by those neighbors within its own block. In terms of polymerization, this could model a situation in which the growing polymers are not evenly distributed throughout the system, but rather, exist in roughly discrete clusters of average size w + 1. The parameters used for both simulations are kBT ) 2.0, fb ) 1, γ0(0) ) 1.0, τ ) 10, ∆t ) 0.007, ∆tξ ) 0.0007, and ζ ) 2. Note that fb is smaller here than in section IIIB, thereby deemphasizing the exoergicity. For the cases with banded-matrix couplings, we set N ) 100; for the block-diagonal cases, N must be evenly divisible by w + 1, so we use N ) (100, 99, 102) for w ) (4, 8, 16), respectively. As a result of the relatively weak biasing force (as compared to that of the previous section), the short-time behavior of the MSVs in the present simulations is qualitatively similar to the free-particle case and hence not shown. The long-time MSVs for all the simulations run with the biased soft-core potential are shown superimposed in Figure 11, indicating that the equipartition result is obtained throughout (except, as noted above, for such short-time deviations as seen in Figure 5). The average position 〈R(t)〉 is shown as a function of time in Figure 12. The equilibrium value of 〈R〉 increases with w, as does the equilibration time. These results are physically appealing when V′sc(R) is viewed as a polymer-growth PMF, since one expects that, all else being equal, longer polymers will interact with more “neighbors” and take longer to form than their shorter counterparts. Note also that the growth profiles appear to be essentially the same for a given value of w, regardless of the specific coupling scheme (banded or block-diagonal) used. The time dependence of the heterogeneity among the local environments, as quantified by the variance in the equal-time value of the friction kernel, is displayed in Figures 13 and 14.
Irreversible Nonequilibrium Environments
J. Phys. Chem. B, Vol. 104, No. 15, 2000 3461
Figure 11. Mean-square velocities for all the simulations using the biased soft-core potential shown overlaid, demonstrating equipartition of kinetic energy obtained throughout. Figure 13. Early-time evolution of the variance in the individual particle environments, 〈γn2(t,t)〉 - 〈γn(t,t)〉2, for the biased soft-core potential at various values of w labeled as in Figure 12.
Figure 12. Average position 〈R(t)〉 displayed as a function of time for the biased soft-core potential at various values of w: (a) w ) 0, (b) w ) 4, (c) w ) 8, and (d) w ) 16. As seen in ref 2, in the w ) N - 1 limit the equilibrium value of 〈R(t)〉 is roughly 8.5, which is lower than that seen in the finite-w cases.
As seen previously for the free-particle and biased-potential cases, at early times the degree of heterogeneity increases with decreasing w, as shown in Figure 13. At longer times, however, the fact that for smaller w the particles are restricted to smaller values of 〈R〉 causes this trend to be reversed in Figure 14. Once again, these observations are consistent with our physical picture of longer polymers which, on average, will sample a wider variety of environments than shorter polymers. Note also that these qualitative observationssas well as quantitative agreement up to fluctuationssare seen between the banded-WiGLE and block-WiGLE results of Figures 13 and 14. This suggests that the size of w, and not the specific composition of the Sw,n, plays a primary role in determining the WiGLE dynamics. IV. Concluding Remarks The self-consistent iGLE has been extended to account for heterogeneity in fluctuating local environments. This results in
Figure 14. Long-time evolution of the variance in γn(t,t) for the biased soft-core potential at various values of w labeled as in Figure 12.
the so-called WiGLEsvis-a-vis iGLE of degree wsin which each tagged coordinate feels a friction due to the dynamics of w other tagged coordinates in its environment. This may be of particular relevance to the discussion of reaction polymerization in which a given polymer may grow in a polycrystallinic or amorphous environment. Notably, in the simulations described in section III C, the equilibrium value of the average coordinate, interpreted here as a mean polymer length, increases with w as we would expect for real polymers: the longer the polymers, the more neighbors, on average, they interact with. In these simulations, we also find that the dynamical effects are dominated by the magnitude of w and not the specific makeup [e.g., banded versus block-diagonal coupling in eq 2.7b] of the neighbor sets Sw,n. In such cases, use of block-diagonal WiGLE can dramatically shorten the computation time by requiring the
3462 J. Phys. Chem. B, Vol. 104, No. 15, 2000 calculation of only N/(w + 1) distinct friction kernels at each time step, rather than the N such calculations required for banded WiGLE. The iGLE, banded WiGLE, and block WiGLE are currently being explored with the polymer PMF of ref 2, in order to better understand polymerization processes. Acknowledgment. This work has been supported by a National Science Foundation CAREER Award under Grant NSF 97-03372. R.H. is a Cottrell Scholar of the Research Corporation, an Alfred P. Sloan Fellow, and Blanchard Assistant Professor of chemistry at Georgia Tech. F.L.S. is a Camille and Henry Dreyfus Faculty Start-up Awardee.
Somer and Hernandez References and Notes (1) (2) (3) (4) (5) (6) (7) (8) 2708. (9) (10) (11) 3172. (12) 1788.
Hernandez, R.; Somer, F. L. J. Phys. Chem. B 1999, 103, 1064. Hernandez, R.; Somer, F. L. J. Phys. Chem. B 1999, 103, 1070. Somer, F. L.; Hernandez, R. J. Phys. Chem. A 1999, 103, 11004. Hernandez, R. J. Chem. Phys. 1999, 110, 7701. Lindenberg, K.; Seshadri, V. Physica A 1981, 109, 483. Carmeli, B.; Nitzan, A. Chem. Phys. Lett. 1983, 102, 517. Corte´s, K. L. E. Physica A 1984, 126, 489. Corte´s, E.; West, B. J.; Lindenberg, K. J. Chem. Phys. 1985, 82, Zwanzig, R.; Bixon, M. Phys. ReV. A 1970, 2, 2005. Metiu, H.; Freed, K. F. Phys. ReV. A 1977, 15, 361. Straub, J. E.; Borkovec, M.; Berne, B. J. J. Chem. Phys. 1985, 83, Straub, J. E.; Borkovec, M.; Berne, B. J. J. Chem. Phys. 1986, 84,
