Molecular dynamics simulation of a polymer chain in solution Burkhard Dünweg and Kurt Kremer Citation: J. Chem. Phys. 99, 6983 (1993); doi: 10.1063/1.465445 View online: http://dx.doi.org/10.1063/1.465445 View Table of Contents: http://jcp.aip.org/resource/1/JCPSA6/v99/i9 Published by the American Institute of Physics. Additional information on J. Chem. Phys. Journal Homepage: http://jcp.aip.org/ Journal Information: http://jcp.aip.org/about/about_the_journal Top downloads: http://jcp.aip.org/features/most_downloaded Information for Authors: http://jcp.aip.org/authors Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions Molecular dynamics simulation of a polymer chain in solution Burkhard DOnwega) [nstitut fiir PhysIk, Johannes-Gutenberg-Universitat Mainz, Postfach 3980, D-55099 Mainz, Germany, and Hochstleistungsrechenzentrum, Forschungszentrum Jiilich, Postfach 1913, D-52425 Jiilich, Germany Kurt Kremer Institutfiir Festkorperforschung, Forschungszentrum Jiilich, Postfach 1913, D-52425 Jiilich, Germany (Received 13 January 1993; accepted 19 July 1993) Results of a molecular dynamics simulation of a single polymer chain in a good solvent are presented. The latter is modeled explicitly as a bath of particles. This system provides a first-principles microscopic test of the hydrodynamic Kirkwood-Zimm theory of the chain's Brownian motion. A 30 monomer chain is studied in 4066 solventparticles as well as 40/4056 and 60/7940 systems. The density was chosen rather high, in order to come close to the ideal situation ofincompressible flow, and to ensure that diffusive momentum transport is much faster than particle motions. In order to cope with the numerical instability of microcanonical algorithms, we generate starting states by a Langevin simulation that includes a coupling to a heat bath, which is switched off for the analysis of the dynamics. The long range of the hydrodynamic interaction induces a large effect of finite box size on the diffusive properties, which is observable for the diffusion constants of both the chain and the solvent particles. The Kirkwood theory of the diffusion constant, as well as the Akcasu et al. theory of the initial decay rate in dynamic light scattering are generalized for the finite box case, replacing the Oseen tensor by the corresponding Ewald sum. In leading order, the finite-size corrections are inversely proportional to the linear box dimensions. With this modification of the theory taken into account, the Kirkwood formula for the diffusion constant is verified. Moreover, the monomer motions exhibit a scaling that is much closer to Zimm than to Rouse exponents (t?-/3 law in the mean square displacement; decay rate of the dynamic structure factor ct:.k3 ). However, the prefactors are not consistent with the theory, indicating that (on the involved short length scales) the dynamics is more complex than the simple hydrodynamic description suggests. I. INTRODUCTION This paper is intended to provide a detailed account on a large-scale molecular dynamics (MD) simulation we have performed on a flexible polymer chain in an explicit bath of solvent particles,1 which is a model system for dilute polymer solutions. Since the late 1970's, this system has continuously attracted the attention of MD researchers,2-8 and has, apart from us, only recently been simulated by two other groups independentIy,9-11 the last paper also providing a nice overview of the historical development. The reason why the single chain in a bath of solvent particles is such a persistent challenge of computational physics is that the problem of polymer dynamics in dilute solution is very amenable to molecular dynamics, but poses a nontrivial computational task: The chain relaxation introduces a long time scale that is much larger than the typical relaxation times of a simple liquid, and that increases strongly with chain length. However, in order to see asymptotic long-chain behavior, rather long chains are required, which, in order to keep the solution dilute, means large system sizes. As will become clear in the sequel, even the present work, which, to our knowledge, is the most extensive computational effort on the problem so far,_has a)Present address: Center for Simulational Physics, Department of Physics and Astronomy, The University of Georgia, Athens, Georgia 30602. apparently not reached the asymptotic limit in every as- pect. As usual, the reason for doing a computer simulation is that analytical theories rely on rather uncontrolled approximations, while experiments are unable to measure all relevant quantities independently with sufficient accuracy. The standard model for the chain's Brownian motion in dilute solution was established by Kirkwood et al. 12,13 and by Zimm,14 who extended the Rouse model15 to include hydrodynamic interactions, which dominate the dynamics in the dilute limit. The equation of motion of this model, Kirkwood's diffusion equation (Ref. 16; also see Sec. II), has been solved analytically, only for the special case of a random walk chain, using the approximation of a preaveraged diffusion tensor. However, in good solvents the chain has the structure of a self-avoiding walk, and, according to Kirkwood's theory, the motion is controlled by a nonpreaveraged diffusion tensor. For this reason, several researchers have used the Brownian Dynamics method to solve the equation numerically.17-22 It should be noted that such an algorithm is rather demanding for long chains, because in every time step the square root ofa 3Nch X3Nch matrix has to be calculated, where Nch is the number of monomers. ___ 9ne can als() exploit the fact that, within the framework of the diffusion equation, there are exact relations between static averages and the short-time dynamics. For J. Chern. Phys. 99 (9), 1 November 1993 0021-9606/93/99(9)/6983/15/$6.00 ®1993 American Institute of Physics 6983 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 6984 B. DOnweg and K. Kremer: Dynamics of a polymer chain instance, the short-time diffusion constant (which is rather close to the long-time diffusion constant), is predicted to be proportional to the inverse hydrodynamic radius (see Sec. II). One can then use a computer simulation to calculate these static averages, and, from that, draw conclusions about the dynamics, taking the Kirkwood diffusion theory for granted. In this case, the method ofchoice is to perform Monte Carlo (MC) simulations,23 which have been very successful in verifying the static scaling laws that govern the chain's conformation statistics.24 This approach of "static dynamics" has been taken, e.g., in Refs. 25 and 26. The data interpretation of dynamic light scattering experiments27-29 usually employs just the opposite approach: While it is very hard to measure, say, the hydrodynamic radius directly, the diffusion constant can be determined from the decay rate of the scattering intensity. From this one can indirectly determine the hydrodynamic radius, taking Kirkwood's diffusion equation again for granted. The present study, however, aims at a test of the validity of the Kirkwood diffusion equation itself. Therefore, Brownian Dynamics or MC cannot be used, because both methods either put the concept into the model from the very beginning, or they do not take hydrodynamic interactions into account at all, because the momentum transport through the solvent is not modeled properly. Strictly microcanonical MD with explicit solvent particles is obviously the most realistic way to do that. As was shown in a previous paper (referred to here as paper 130), even an MD with noise results in unrealistic dynamics. An alternative approach to MD are lattice-gas cellular automata (LGCA),31 or hybrid schemes between LGCA and MD.32 These are modern and very powerful techniques, which, however, are not yet fully understood. Mainly due to limitations in computer resources, the older MD studies2-8 have not been able to accomplish the goal set in the previous paragraph. These simulations were confined to rather short chains, and mainly static quantities were analyzed, for reasons of statistical accuracy. The present study was done with a similar computational effort as a previous large-scale MD simulation on the dynamics of polymer melts.33,34 This enabled us to quantitatively compare the dynamical properties of the chain with the predictions resulting from the Kirkwood diffusion equation. In summary, we find very good agreement of theory with computer experiment on long length scales, the Kirkwood prediction for the diffusion constant being very nicely verified. On shorter length scales, we find that the scaling predictions ofthe Zimm model (r/3 behavior ofthe mean square displacement, k3 t decay of the dynamic structure factor) hold for our system. However, the prefactor is not consistent with the prediction resulting from the hydrodynamic theory. Therefore we argue that these shOIter length scales are already too close to atomic length scales, and the simple hydrodynamic picture breaks down. The paper is organized as follows: In Sec. II, we briefly summarize the theoretical predictions of the Zimm model and introduce the notation. Section III contains a detailed discussion of the expected finite system size effects, based on the theory of Ewald sums. In particular, we present an analytical calculation (with severe approximations), which yields a closed expression for the box size dependence of both the diffusion constant and the initial decay rate of the dynamic structure factor. In Sec. IV we define the simulated model and describe the simulation procedure we applied. In Sec. V we present the results for the solvent properties, with emphasis on those that are important for the chain dynamics. In Sec. VI we discuss the most important static properties of the chain, while Secs. VII and VIII contain the results for its dynamic properties (diffusion constant and local motions, respectively). In Sec. IX we conclude with a discussion. II. THEORY FOR THE IDEAL SYSTEM Consider a single long flexible polymer chain of Nch monomers in an infinite bath ofsolvent particles in thermal equilibrium. The solvent shall be good; hence the static configurations are characterized by the exponent v~O.59 that relates lengths in real space with lengths along the chain, e.g., the end-end distance, (R2) = «rNch -rt)2), and the radius of gyration, (R~)= 2N\ ~ (";j) ch IJ (1) (2) (rij=rj-rj' rj being the monomer coordinates), scale with the chain length as (R2) ex: (R~) ex:N~~. (3) Similarly, the static structure factor, -I " -I " (Sin(krij )) S(k)=Nch ~ (expUk.ri)=Nch ~ kr·· ' IJ IJ IJ_ (4) which is measured in scattering experiments, in the regime Rat 00 and ka .....O. For a real system, the finiteness ofRG and a both tend to decrease D(k) compared to Eq. (23). A more detailed consideration of these effects shall be published elsewhere.38 The dynamic scaling relations are summarized as follows: In the Rouse case (which we consider too, in order to illustrate the influence ofhydrodynamic interactions on the dynamic exponents), one has D N -I R-lIva: ch a: G , while with hydrodynamic interactions D R-I R-I a:Ha:G' (24) (25) It should be noted that the hydrodynamic radius has very large corrections to scaling,26 which can be understood in terms of finite chain length and finite bead size.38 The longest relaxation time l' (which is called Rouse time 1'R or Zimm time 1'z, respectively) can be estimated via 7=R~(6D), (26) which is the time the chain needs to move its own size. Hence, (27) and 1'za:Rb· . _.- (28) This defines the dynamic exponent z=2+ l/v:::::3.7 in the Rouse case, while it is z=3 in the Zimm case. Note that without hydrodynamic interactions, z depends on v (I.e., the dynamic scaling is different for chains in good solvents and in e solvents), while with hydrodynamic interactions included, it is independent of the solvent quality. This exponent also appears in the subdiffusive behavior of the mean square displacement of a monomer on time scales intermediate between microscopic times and the longest relaxation time: (29) The exponent is :::::0.54 in the Rouse case ·and -~ in the Zimm case. Similarly, the dynamic structure factor obeys the relation (on the same time scales and for R(jl<,k<,a- I ) (30) III. THEORY OF FINITE SYSTEM SIZE EFFECTS Due to the long-range nature of the hydrodynamic interactIon one must expect strong effects on the dynamical properties by the finiteness of the MD box, even if the chain fits very nicely into the box and the static properties are not affected at all. Since the simulation is run with periodic boundary conditions, the chain has an infinite number of periodic images with which it interacts as well. ·One might instead consider running the simulation with, e.g., hard walls, but such a modification would not remove the finite-size effect, but just replace it by another effect, which is much less controlled and understood, and probably much larger. There is no other way of removing the finite-size effect than a sufficiently large system, which is computationally not feasible. A good and quantitative understanding of the influence of the image chains is therefore essential. Intuitively, the effect can be viewed as an effective increase of the hydrodynamic radius due to the image chains causing the diffusion constant to decrease. It should be pointed out that this is a completely different mechanism than an increased concentration of the dilute solution: In a dilute solution of many chains the relative positions of the chains vary strongly, while they are fixed in space for the present situation. To put it more formally: In the case of the periodic box, the number ofdegrees offreedom appearing in the Kirkwood diffusion equation is not increased by . theadditional chains, while in a multichain system it is. Quantitatively, the effect can be taken into account by modifying the diffusion tensor appropriately. Denoting the linear size of the cubic MD box with L, one has for i=l=j (cf. the previous paper30), Oij=O(rij), with kBT " l-klSlk OCr) =~L ~ ,.2 exp(zk'r), 'TJ k#O K- (31) (32) where k= (21TU)/L runs over the reciprocal lattice vectors of the MD box [n is a vector ofintegers). This form comes from the superposition of hydrodynamic modes. From the k-2 dependence one sees that the main contribution is due to the long-wavelength modes. k=O is excluded because this would correspond to an overall translational motion of the system, and we study the diffusion of the chain relative to the fluid, which is globally at rest. While in the infinite system all modes down to k=O contribute (the sum is replaced by an integral), in a finite box the longwavelength modes are cut off, resulting in a slowing down of the diffusivity. Similar finite-size effects have also been found in MD simulations of solids, where the finite system size cuts off the long-wavelength part of the phonon density of states.39 Note that the explicit exclusion ofk=O in Eq. (32) is necessary, not only for the physical reasons mentioned above, but also mathematically, in order to keep the expression finite. The series converges because of the oscillations of the exponential. An attempt to write down the J. Chem. Phys., Vol. 99, No.9, 1 November 1993 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions B. DOnweg and K. Kremer: Dynamics of a polymer chain 6987 hydrodynamic interaction with the periodic images by just adding up the infinite number of Oseen tensors in real space would lead to a mathematically ill-defined formula, because this series contains a k=O component and hence diverges. For this reason, the Fourier representation, Eq. (32), is more useful. The diagonal elements also have to be modified due to the hydrodynamic interaction of a bead with its own periodic images. This can be done via Dii=Dol+lim(D(r) _skBT (l+P® P»), (33) r ...O 1T1V where D (r) according to Eq. (32) is used. For numerical purposes, the formulas, Eqs. (3Z) and (33), are not very useful because the convergence of the series is very slow. However, the sums can be rewritten in a rapidly converging form, using the Ewald summation technique.40 For the Oseen tensor, formulas have been given by Beenakker.41 In the Ewald representation it is also easy to perform the limit r->O in Eq. (33). Numerical evaluation of the remaining series then yields 1 2.S37kBT 3Tr Dii=Do 61TTJL' (34) where the constant -Z.837 is the analog of a Made1ung constant. This allows us to define an effective, box-sizedependent hydrodynamic radius via the relation 2.837kBT 1 " 6 LN +3N2 .~. Tr(Dij) 1TTJ ch ch z=/=J Do kBT / 1 ) =Nch+61T1J \RH L' (35) Di} again taken from Eq. (3Z) to include the periodic im- ages. From this result one explicitly sees that the chain diffusion constant is affected by the finite-size effect, and that therefore the Zimm time is changed. In a simple-minded finite-size scaling theory that we presented in Ref. 1, we assumed that the effect on the other decay rates of S(k,t) is just a change by a uniform factor. Pierleoni and Ryck- aertlO tested this assumption by MD simulation of short chains in boxes of various size and found it not supported by the data. That it is indeed wrong can also be seen analytically by inserting the finite-size corrected diffusion tensor in Eq. (Z2). Although an exact evaluation is only possible numerically, a simplified and approximate treatment that shows the main features can still be done. First, we rewrite Eq. (22) in the same spirit as Eq. (4.110) of Ref. 16 as k T i l J 1 (k' QA)2 D(k) = TJ~3 S(k) 41T d 2 k Q~ (f S(k-Q). (36) Here the Fourier representation of the Oseen tensor, Eq. (32), and the definition of the static structure factor, Eq. (4), have been used. The finite system size shows up by the discrete sum over the wave numbers Q= (21711)1L. Moreover, an explicit spherical average [represented by (41T)-1 Jd2 k] has been introduced in order toremove the cubic anisotropy. Using the spherical symmetry of S(k) (the statics is unaffected by the finite box size) and introducing u= k.Q, this is rewritten as kBT I·S·1 D(k)=-::y; L -2· du(1-u2) TJ Q*O -I S( ~~+Q2-ZkQu) X (fS(k) (37) In order to make this easily tractable, we replace the sum over Q by a spherically symmetric integral, and take the finite-size effect approximately into account by restricting the integration volume on the whole Q space, except a sphere around Q=O with radius (21Tao)/L, where ao is a parameter of order unity. As a further simplification, we consider the structure factor for a random walk chain in the approximate form16 S(k) Nch (38) The Q summation is then very easily done with the result 1 kBT fl 2 I+K2/Z D(k) =Z-2 -R du(1-u ) ~ 2 2 1T TJ G -I Z+K (l-u ) _X (i+arctan ~2~~~~U2»). (39) Here the abbreviations K =kRG and E= (Z1TaoRG)1L have been introduced. The chain diffusion constant is calculated from that for the special case k=O, with the result D=311T~~ [1-~arctan( ~1Tao:G)], (40) showing that the finite size effect decreases the diffusion constant, controlled by the parameter RoiL. In the scaling regime K> 1, Eq. (39) is simplified to I kBT flD(k)=4rTJR G K -I du~l-u2 X [i+arctan(U~~~~)]. (41) Since KIRG=k and EIK= (Z1Tao)l(kL), one sees that the dependence on RG drops out. Instead ofRoiL, the relevant dimensionless parameter for the finite-size effect is now 1!(kL). Linearizing the above result [Eq. (41)] with respect to E, one obtains, for large L, 1 kBT 2ao kBT 2 D(k)=16nk-31T TJL +O(L-). (42) The first term is identical to the well-known result for D(k) for a random walk chain in the asymptotic limit kRG>1, ka ~ 1.16 ,36 The second term gives the finite-size effect in the scaling regime, which is here predicted as a constant, k-independent shift. This shift is the same as J. Chern. Phys., Vol. 99, No.9, 1 November 1993 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 6988 B. D(jnweg and K. Kremer: Dynamics of a polymer chain what one obtains from Eq. (40) upon linearization with respect to L -1, i.e., in linear order of inverse system size the shift is the same for all k values. This can also be shown directly from Eq. (39) by linearizing with respect to the expansion parameter E. In principle, the results Eqs. (40) and (42) mean that the finite-size effect corrupts the scaling behavior both in the kRG~ 1 as well as in the kRd~ 1 regime: The additive shifts cause the laws Da::.R(;I and D(k) a::.k, respectively, to be no longer strictly valid. Moreover, the result shows that the relative contribution of the finite-size effect gets weaker when k is increased. This is so because D(k) increases with k, and qualitatively consistent with the observations of Pierleoni and Ryckaert. 1O This, however, shows that the way how they interpreted their data is not the only possible explanation why finite-size effects become less and less pronounced with increasing k. They suggested a completely different mechanism based on a retardation argument: The finite-size effect simply cannot build up on the time scale of the decay of the structure factor, because the hydrodynamic interaction needs a finite amount of time to spread to the periodic images. This effect should become more and more pronounced with increasing k, because there the structure factor decays more rapidly, and because more time for the spreading is needed, since "domains" of size 21T/kare further apart from their images. This argument is intuitively rather appealing, however, it is not completely obvious to us that it is indeed correct. A more rigorous consideration based on a retarded interaction tensor in the framework of Kirkwood's diffusion equation is certainly desirable. Here we only want to point out that the empirical observation of a decreased finite-size effect at higher k does not prove the retardation argument, since, according to the above calculation, one has to expect the same behavior, even in the limit of no retardation. Quantitatively both their and our MD data for D(k) at the higher k values (cf. Sec. VIII) disagree with the theoretical predictions-regardless of whether one uses the theory for L -> 00 (this would be correct if the PierleoniRyckaert retardation mechanism suppresses the finite-size effect) or the finite-size corrected predictions. Therefore we assume that both their and our data are outside the hydrodynamic regime, as we will discuss in more detail in Sec. VIII. IV. THE MODEL AND SIMULATION TECHNIQUE The first step of constructing a MD model for the dynamics of dilute polymer solutions is the definition of the solvent in terms of explicit particles. As has become clear from the discussion of the theory, the solvent mainly exists to transport momentum. Therefore, we need particles with a short-range strong repulsive interaction. An attractive tail may be included, but is not necessary. Therefore we only used purely repulsive interactions in order to save computer time. The "ideal" solvent would be a dense hardsphere liquid. However, the project was performed on a vector computer (Cray-YMP at the HLRZ Jiilich), and for hard spheres, there is no known fast vectorizing MD algorithm. For continuous potentials, on the other hand, the program can be efficiently vectorized using the layered link-cell method,42 which we applied here. Hence we used "soft spheres" interacting via a truncated Lennard-Jones (WCA43) potential: Uu(r) = 14e [ (~r2-(~r+~], r<2 116 u, (43) 0, r>2116 u. Together with the particle mass m, this potential fixes the unit system: Lengths are measured in units of u, times in units ofru=(mu2/e)1I2, maSSes in units ofm, energies in units of e, etc. To simplify notation, we will henceforth mostly use a unit system in which all three parameters are unity. We simulated the system at temperature kBT= 1.2 and at density p=1.05-3 ::::;0.864. This is a very high density, which was chosen intentionally in an attempt to match the ideal fluid of the theory most closely: The theory assumes incompressible flow, and one should expect a low compressibility at high densities. In order to analyze the solvent properties, we, at first, ran the pure solvent. The particles were set up on a simple cubic lattice (which quickly melted) and equilibrated either by assigning initial random velocities, followed by velocity rescaling to fix the kinetic energy, or by coupling the system to a viscous background via friction and random force,30,44 until the original lattice correlations completely had decayed. Dynamical properties were analyzed, after switching off the equilibration procedure, in purely microcanonical runs. In view of the analytical result of paper 130 on the modification of the time· correlation functions by noise, we regarded any deviation from strictly Newtonian dynamics as too dangerous. Typically, we used a time step of 0.004 in the runs without heat bath, and 0.006 with heat bath. In cases where we wanted an accurate estimate for the short-time behavior (e.g., in order to measure the viscosity), we used a time step of 0.001. A fifth-order predictor-corrector scheme45 was applied to numerically integrate the equations of motion. Without heat bath, at a time step of 0.004, our program42 performed 3X 105particle updates per second on a single Cray-YMP processor, independent of the system size, which varied in the range from 93= 729 to 203= 8()()Q particles. A polymer chain was introduced into the system by redefining the first Nch out of the N tot particles as monomers that were connected by an attractive backbone potential. As an initial configuration of the chain, we generated a self-avoiding random walk on the initial simple cubic lattice of particles. For the backbone potential we added the FENE potential, k 2 ( ?)Uch (r)=-2"Roln l-R~ , (44) with the parameters k=7 and Ro=2 to the repulsive interaction, Eq. (43). The nonlinearity restricts the bond length to Ro at maximum. No bond angle dependence was introduced in order to make the chain as flexible as possible. Note that we made the potential much "softer" than J. Chern. Phys., Vol. 99, No.9, 1 November 1993 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions B. DGnweg and K. Kremer: Dynamics of a polymer chain 6989 the corresponding potential in the melt studies.33,34 We did this in an attempt to increase the coupling of the monomers to the surrounding solvent relative to the coupling along the chain. Similarly, we increased the monomer mass to the value 2 in order to make them behave more like Brownian particles. The effect of mass change on the diffusion of a particle in a surrounding of otherwise identical particles has been studied by Toxvaerd,46 using MD. However, it is not entirely clear how far these modifications changed the dynamical properties of the chain, particularly on long length scales. Although these are very interesting questions, we did not try to systematically study them, because this would have required runs for different models beyond the scope and CPU time of the present investiga- tion. Otherwise there is no difference between solvent particles and chain monomers. In particular, the repulsive Lennard-Jones interaction acts in exactly the same way between all particles. Apart from simplifying the simulation program, this feature also removes uncertainties about the theta transition: The solvent is ideally good (a so-called "athermal" solvent24 ), and the theta collapse never occurs. Therefore, the good-solvent condition is the natural choice that is most easily modeled in a computer simulation. We studied three systems of a chain in solvent, an Nch=30 chain in Nso1v=4066 solvent particles, as well as Nch=4OINsolv=4056 and Nch=60INsolv=794O systems. In order to explore the configuration space of the chains, we ran them for 2.1X105 'TLJ (Nch=30), 2.238X105 'TLJ (Nch=4O), and 1.95 X 105 'TLJ (Nch=60). This means more than 30 million integration steps, using a time step of 0.006. On this time scale, we found ourselves unable to keep a purely microcanonical MD numerically stable (the instability is due to the dicretization errors caused by the finiteness of the time step). Therefore, we coupled all particles to a heat bath by a MD with noise, using a friction constant of ~=0.5. As discussed in paper 1,30 this method generates a canonical ensemble, but is in serious error as far as the hydrodynamic properties are concerned. We therefore used these runs only in order to generate initial states, saving the system configuration every 300'TLJ. Mter that, dynamical properties (i.e., time correlation functions) were obtained by averaging over (roughly 7(0) runs without noise, which started from these configurations, and then ran for l00'TLJ with a time step of 0.004. This dynamical observation time is significantly shorter than the Zimm time (for Nch=30, it is roughly O.4'Tz, while for Nch=60 it covers only about O.I'Tz). However, it is long enough to reach well into the interesting intermediate time regime. The first 20% of these runs were discarded in the analysis, in order to allow for the building up of the hydrodynamic correlations beyond the screening length I of the equilibration run. According to paper 1,30 the value of I is 1=&~2.4' (45) where we used the viscosity 7]=2.4 (Sec. V). On the other hand, the radii of gyration are 3.28, 3.78, and 4.78 for 3.0 r----.----------,----, 2.0 1.0 0.0 ~~-...LL.-~---'---~--'----.J 0.0 1.0 2.0 3.0 FIG. 1. The pair correlation function of the liquid as defined in Eq. (46). Nch=30, 40, and 60, respectively. Hence one expects a strong screening effect, in particular, when the chain length increases. Indeed, a rough estimate of the chain diffusion constant in the equilibration runs yields D=4.8X 10-3 , 3.4X 10-3, and 2X 10-3 for Nch=30, 40, and 60, while without friction we measured the corresponding values D=6.82X 10-3, 5.45 X 10-3, and 4.25 X 10-3. In order to assess the statistical quality of our procedure, we estimate the longest relaxation time in the equilibration run via Eq. (26), yielding 1'=370, 700, and 1900 for Nch=30, 40, and 60, respectively. (Note that this is not 'Tz!) The observation times of the canonical runs were larger by a factor of 570, 320, and 100, respectively. V. SOLVENT PROPERTIES The pair correlation function of the liquid for r> 0 is defined as g(r) (p(r)p(O) ) (p)2 where per) is the local density, per) =L-3 L8(r-ri)' i (46) (47) Because of spherical symmetry, g depends only on the distance r. For our system, g(r), which is shown in Fig. 1, exhibits a very pronounced first-neighbor peak and welldefined second- and third-neighbor shells, indicating that our solvent is highly correlated and exhibits structure at least up to length scales of r~ 3. This is not surprising, in view of the high density, which shows also up in the high value of the pressure P=9.84 with only little fluctuations, ~(ji2) - {p)2 = 0.03. Pwas evaluated as P=~Tr P, (48) where the pressure tensor P is obtained via the virial the- orem:47 .(49) J. Chern. Phys., Vol. 99, No.9, 1 November 1993 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 6990 B. DOnweg and K. Kremer: Dynamics of a polymer chain 3.0 ,--~-____r----,----_.__--___, 2.0 1.0 0.0 '---~--'---~-"'---~-"'---~-----' 0.0 0.5 1.0 1.5 2.0 FIG. 2. 7](t) for an 8000 particle system. where Fij stands for the interparticle force between particles i and j, and a, /3 denote Cartesian indices. The offdiagonal elements of P can be used to determine the shear viscosity 71 via Green-Kubo integration:47 71= L3 . roo dt(p43(0)pf3(t), kBT Jo . (50) for a=l=/3. We performed the integration for various system sizes (up to 8000 particles), using runs with time step 0.001 of duration 1OOOru, and obtained the estimate 7J=2.4±0.1. No finite-size effect was observable within the statistical error. The function 7J(t), which is just the integral up to time t, is displayed in Fig. 2 for the largest system. The kinematic viscosity, which governs the diffusive momentum transport, then results 7Jldn=2.8±0.1. From the same runs, we also studied the diffusion constant Dso1v of the particles via their mean square displacement. Since the simulation included a very slight drift of the overall system (e.g., the 8000 particle system had a drift velocity of 2.67 X 10-2 ), we calculated the mean square displacement relative to the center of mass. Moreover, the different systems did not run at precisely the same temperature (as measured by the kinetic energy): Typically, the actual temperature differed from the "ideal" value 1.2 by about 0.5%. In order to assess the (rather small) finite-size effect, we therefore had to study the mobility J-L=Dso1vl(kBT). According to Eq. (34), we should see a L -I behavior with slope -2.837/(61T7J). Figure 3 shows our result. There is some statistical scatter in the data, but qualitatively the L -\ behavior is well confirmed. The line is a regression fit with slope -0.0593, which corresponds to an effective viscosity of 2.54. Within our errors, we may say that the prediction of Eq. (34) is confirmed by the data. From the regression fit, we also obtain the extrapolated diffusion constant for the infinite system, Dso1v =0.078, i.e., the kinematic viscosity is indeed much larger than the particle diffusion constant, meeting the consistency requirement, Eq. (15) by a factor of 36. This means that the transport properties are strongly dominated by momentum transport instead of mass transport. More- 0.063 ,---...-'---,---_--'-....,...:.-_---,--_----, 0.062 • O.06j •;; 0.060' 0.059 . -- 0.058. ~~~--::::::_-~-=-=-~-...,....,,::__...,....,.--=-'. 0.04 0.06 0.08 0.10 0.12 1/L FIG. 3. The solvent particle mobility /L as a function of inverse linear system size L -I. The straight line is a regression fit to the data. over, we note that the prediction from a Stokes law with slip boundary conditions, Dso1v = (kBT)/(41T7Ja) , yields a value Dso1v=0.075 if we use 71=2.4 and a=0.53 [the latter estimate as a typical particle radius, since the typical nearest-neighbor distance, estimated from the first maximum of g(r), i81.06]. Finally, we also calculated the time autocorrelation function of transversal velocity field modes UkA., as defined in paper r,30 (51) where N is the number of particles, Pi denotes the particle momenta, k is the wave number of the mode, and EA. is a polarization unit vector orthogonal to k. Since this calculation is rather costly, we did it only for one system of2197 particles, choosing k= (21Tnez )/L, n= 1,...,10. The hydrodynamic theory predicts (cf. paper r30) -- roo ~T D(k):= Jo dt(ut.(O)UkA.(t)=~. (52) This is of direct interest, since this relation enters the derivation of the hydrodynamic interaction; essentially jj(k) is just the Fourier transform of the Oseen tensor. For long wavelengths the relation should hold, while for higher k values one has to expect deviations. This gives a clue about the length scales down to which the hydrodynamic Oseen description is expected to hold. Figure 4 shows our result: For the long-wavelength modes, the ~reement is indeed excellent, while for the higher modes D(k) is significantly above the hydrodynamic prediction. Here the atomic length scales come into play, and the autocorrelation function of UkA. decays more slowly. For later, we note the following important implication of the result: If one attempts to remedy the deficiencies of the Oseen tensor on short length scales by introducing a generalized diffusion tensor, which, for long distances, asymptotically tends to the Oseen tensor, then the most I!~.tural generalization would be the Fourier transform of D(k), as measured in the computer experiment. This would still assume (cf. paJ. Chern. Phys., Vol. 99, No.9, 1 November 1993 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions B. DOnweg and K. Kremer: Dynamics of a polymer chain 6991 10" g10-4 Cl • •••••• 10·' 10° 10' k FIG. 4. Log-log plot of D(k) [Eq. (52)] according to the Green-Kubo definition (points), and according to the hydrodynamic prediction, as stated by Eq. (52) (straight line), using a viscosity value of 1/=2.4. Note that there is no free parameter to adjust the theory to the data. per 130 ) that the monomers are dynamically "embedded" in the solvent surroundings (so that one can just use the solvent correlation function in order to calculate the Green-Kubo integral for the monomers), and that only transversal modes contribute. As we will see later, our polymer data indicate that such an attempt is bound to fail. VI. STATICS OF THE CHAIN For reference, we list in Table I some important static and dynamic data characterizing the chains. We measured the gyration radii RG= (R~) 112 as 3.28, 3.78, and 4.78 for Nch=30, 40, and 60, respectively. The corresponding box sizes are 16.8, 16.8, and 21, indicating that the chain fits well into the box and that self-overlap effects are small. Hence we should see unperturbed goodsolvent behavior. Fitting a power law, RG=A (Nch- 1)v, (53) to the data yields A =0.54 and v=0.53. Similarly, we measured the end-end distances R=(R2)1I2 (8.39,9.41, and 12.35 for Nch=30, 40, and 60, respectively), and applied TABLE 1. Summary of some chain properties: Number of monomers Nch ' linear system size L, radius of gyration R G , end-end distance R, inverse hydrodynamic radius Rill (both for an infinite solvent as well as for the finite MD box), diffusion constant D, Zimm time 7Z, and monomeric diffusion coefficient Do, as obtained from comparison with the Kirkwood formula for D. Nch 30 40 60 L 16.8 16.8 21.0 RG=(R~}112 3.28 3.78 4.78 R=(R2) 112 8.39 9.41 12.35 RaiL 0.195 0.225 0.228 (Rill) L~", 0.34 0.30 0.25 (Rillh 0.18 0.15 0.13 D 6.82XlO-3 5.45 X 10-3 4.2SxlO-3 7z::::R~(6D) 260 440 900 Do 0.062 0.062 0.055 10° k FIG. s. Log-log plot of the static structure factor of the chain, for Nch =30 (lower curve), 40 (middle curve), and 60 (upper curve). the same analysis with the result A =1.3 and v=0.55. This indicates that for both radii the asymptotic long-chain limit is approached. A much clearer picture is obtained from the static structure factor, Eq. (4), displayed in Fig. 5, which we evaluated using the form on the right-hand side for an optimal spherical average. The k- IIv decay is already well established for a sufficiently large window of wave numbers, and can be used to determine an effective exponent v=0.58±0.01. For the inverse hydrodynamic radius according to Eq. (17) we obtained the values (RIll) 00 =0.34,0.30, and 0.25 for Nch=30, 40, and 60, respectively. The corresponding finite-size corrected values (RIll) L according to Eq. (35) are 0.18, 0.15, and 0.13, illustrating the strong influence of the finite box size. The dependence of the effective hydrodynamic radius on the box size is demonstrated in Fig. 6 for the Nch=30 chain, where (RIll) L is plotted versus L -1. These numbers were obtained by putting the chain configurations into "virtual" boxes of various sizes, and 0.4 ,-------.-~---,--_---r--~-, ~..0.3 • •1\ "":":r: 0.2 a: • /' v 0.1 • •0.0 l..--~_--'-_~_-'-_~_-'--_~___It 0.00 0.05 0.10 0.150.20 L·1 FIG. 6. The effective inverse hydrodynamic radius (Rill> L as a function ofinverse linear box size L -I, for the configurations ofthe Nch = 30 chain. The condition of the simulation (L=16.8) is indicated by an arrow. J. Chem. Phys., Vol. 99, No.9, 1 November 1993 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 6992 B. DOnweg and K. Kremer: Dynamics of a polymer chain 0.08 0.07 •0.06 • •0.05 :::J" ~ 0.04 •£:) 0.03 0.02 0.Q1 0.00 0.0 0.1 0.2 0.3 0.4 0.5 L·1 FIG. 7. The finite size effect of D(k), as defined in Eq. (21), and evaluated with the static average, Eq. (22), using Ewald sums. We used the configurations of the Nch=30 chain at k=3.75. The system size of the simulation (L= 16.8) is again indicated by an arrow. doing the Ewald sum analysis for them. One sees that the exact evaluation of the Ewald sum yields qualitatively the same functional form as the approximate calculation, Eq. (40). Moreover, the figure also shows that a box that would be large enough to reduce the finite-size effect on RH down to 10% would have to contain of the order of 3X 105 particles! In the same spirit, we also used the static configurations in order to calculate the Akcasu et al formula for the initial decay rate of the dynamic structure factor, Eq. (22). The denominator is (except for normalization) just the static structure factor, and was evaluated, as described previously. In the numerator, we inserted the Ewald sum of the Oseen tensor for Dij , using the parameters Do=0.06 (see the next section) and 71=2.4. Since the finite system size introduces a cubic anisotropy, the comparison with the dynamics data becomes more meaningful if the numerator is explicitly spherically averaged [S(k,t) was obtained using the explicitly isotropic formula on the right-hand side of Eq. (20)]. Therefore we generated 30random k directions uniformly distributed on the unit sphere, and used them for an approximate spherical average of the numerator. This evaluation is relatively expensive, and hence only a limited number of parameters was used. Analogously to Fig. 6, Fig. 7 demonstrates the finitesize effect on D(k) for k=3.75. The Ewald sum was evaluated for different sizes of the virtual box into which we put the configurations of the Nch = 30 chain. As anticipated by the approximate analytical calculation of Sec. III, one sees that at this higher k value the relative contribution of the finite-size effect is much smaller than at k=O. Figure 8 shows the Ewald summation result for D(k) as a function ofk, where we used the actual box sizes of the simulation runs. The data for the three chain lengths are rather close to each other, indicating that in the studied range of k the effects of both finite chain length and finite system size are not very important. However, the data substantially differ from the asymptotic Akcasu-Benmouna ··0.06 2" 0 0.04 0.02 - 0.00 ~~--L_~-'-_ _-'-~_-,--~_,--~--, o 2 3 k 4 5 6 FIG. 8. D(k) [Eq. (21)] according to the Ewald sum evaluation ofEq. (22), for Nch =30 (points), Nch=40 (triangles), and Nch =60 (diamonds). The straight line is the asymptotic Benmouna-Akcasu formula, Eq. (23), using a viscosity value of 1/=2.4. formula, Eq. (23), which is also included in the figure. We think that this is mainly an effect of the finite size a of the monomers; Eq. (23) is only valid in the limit ka.....O and kRG..... 00. Indeed, a simple model calculation,38 which takes the finiteness of both RG and a into account predicts shifts due to both nonidealities, which both tend to decrease D(k). VII. THE CHAIN DIFFUSION CONSTANT Figure 9 shows the behavior of the mean square displacement of the chain's center of mass as a function of time t. The diffusion constant D was obtained from this by fitting a straight line in the time interval [10,60], yielding the values D=6.82X 10-3, 5.45 X 10-3, and 4.25 X 10-3 for Nch = 30, 40, and 60, respectively. Times t < 10 were not taken into account, since on these time scales the ballistic behavior had not yet died out completely. The data t> 60 were disregarded for reasons of statistical accuracy. 4.0 3.0 A '"c: 2.0 giving a value of about 0.06 for all three systems, which is about 20% smaller than the particle diffusion constant of the solvent. VIII. RESULTS ABOUT MONOMER MOTIONS In Fig. 10 the time dependence or- the mean square displacement of the single monomer is shown in a log-log plot. «ar)2) was averaged over only the four innermost monomers, in order to minimize end effects. A comparative evaluation of four monomers at the chain end indeed indicated a slightly steeper slope of «ar)2) vs t. Similar "'~ 10° 0.05 (above noise level), and to the scaling regimes 20<:;;t<:;;80 and 0.7<:;;k<:;;3, as obtained independently from «ar)2) and S(k), respectively. In Figs. 12 and 13, we compare the Rouse and Zimm predictions with each ""'~ O there is nice agreement between dynamics data and theoretical prediction (for k=O this is just the result of the previous section, and for k=O.5 it is visible in the figure). However, for the higher k values the actual - 0.08 r--~----,------r-----_.__-___; 0.06 z-o 0.04 0.02 2 4 k s • 6 FIG. 15. D(k) as defined in Eq. (21), obtained from the dynamical data for Nch=30 (filled circles), 40 (filled triangles), and 60 (filled diamonds). Instead oftrying to perform the limit t....Owe took the maximum value 9f D(k,t) [cf. Eq. (21)]. For comparison, the data resulting from the static evaluation with Ewald sums (cf. Fig. 8) are also included with corresponding open symbols. dynamics is slower_than the Akcasu-Benmouna predic- tion. In order to try to understand this, let us recall the basic assumptions that led to the Akcasu formula, Eq. (22) (also cr: paper 130). (1) The dynamics is described via a diffusion equation on sufficiently long time scales. (2) The monomers are assumed to be "embedded" into the solvent flow field, and hence their velocity autocorrelation function is replaced by that of the solvent. (3) Only transversal hydrodynamic modes are taken into account. (4) In a long-wavelength approximation, D(k) [cf. Eq. (52)] behaves as k-2 , corresponding to the r-1 behavior of the Oseen tensor. .Assumption 1 is probably well justified. Moreover, we think that the reason for the discrepancy cannot be a breakdown of assuption 4 alone: If this were the case, one could correct for that and instead insert a generalized Oseen tensor, where for D(k) the long-wavelength approximation is not used. However, according to our result, Fig. 4, a more realistic choice of the interaction tensor would increase D(k) compared to the pure Oseen behavior, and henceD(k) according to the Akcasu formula, Eq. (22), would be increased, too. Such a correction would therefore cause an even larger discrepancy. We therefore think that on the length scales in consideration, assumptions 2 and 3 must be questioned as well. The discrepancies seem to become important at k-z 1, corresponding to a length scale of r= (21T)lk-z6. This can no longer be viewed as large compared to atomic length scales: g(rl (Fig. 1) exhibits some structure even beyond r=3, and D(k) (Fig. 4) shows deviations from the purely hydrodynamic behavior for k> 1 as well. It seems quite conceivable that on these length scales the chain can no longer be considered as embedded into the surroundings, but dynamically starts to "emanciJ. Chem. Phys., Vol. 99, No.9, 1 November 1993 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions 6996 B. DOnweg and K. Kremer: Dynamics of a polymer chain pate" from the solvent, with which, of course, it still remains strongly coupled. Pierleoni and RyckaertlO interpreted their data for the decay ofS(k,t) at short times and high k values as asymptotic Zimm behavior, which could be used to extrapolate to longer chains. As discussed at the end of Sec. III, they suggested a retardation mechanism that would suppress finite-size effects. We have severe doubts concerning this interpretation. The problem is that the hydrodynamic interaction does not only need a finite amount of time to spread over the system, but it also needs time to build up, i.e., "forget" the ballistic short-time behavior. Our particles needed of the order of 10-20Tu, until they reached the Brownian diffusive limit. Therefore, we think it is questionable if it is justified to regard the shdrt-time regime of S(k,t) down to O.5Tu as the correct regime to compare with the theory (even if one takes into account that the simulation of Ref. 10 was run at a different liquid state point). Moreover, these authors see simila~ discrepancies in the prefactor of D(k) as we do.49 Altogether, it seems that the scaling laws D(k) ex: k and S(k,t) =S(k,O)f(12t) hold down to both length and time scales, where there is no longer any known theoretical reason to assume that they would have to. However, it has also become evident that the validity of the scaling down to these scales does not mean the validity of the KirkwoodZimm model down to these scales-otherwise the prefactor of D(k) would have to coincide with the theoretical predictions. We think there must rather be another mechanism that causes the behavior. This might well be a complicated interplay between hydrodynamic effects and atomic motions; at any rate it is something "beyond Oseen." We think this regime is, in essence, not understood and poses a very interesting challenge for transport theory. IX. SUMMARY AND OUTLOOK The present work demonstrates that today MD simulations are able to successfully attack the problem of Brownian motion in complex fluids from a first-principles point of view without a priori assumptions about the hydrodynamic interaction. It shows viable routes how to cope with the technical problems of long relaxation times and the large influence of the finite system. size on the dynamical properties. Suspensions and semidilute polymer solutions can now be treated similarly. However, not even the simple case of dilute polymer solutions is fully understood yet on the small length scales of the simulation: While we find very good agreement between dynamical data and hydrodynamic prediction on long length scales [D(k) for k->O], the value of D(k) deviates considerably from the theoretical prediction for the higher k values, indicating that for atomic length scales the dynamics is more complex. Nevertheless, scaling still seems to hold in this regime. A theory providing a description and explanation of these observations would be extremely valuable, in particular also for the interpretation offuture MD simulations of polymer solutions. We feel that a test/reproduction of every aspect of the Zimm model by MD simulation in the fully asymptotic regime, where the chain size is much larger than the correlation length of liquid structure, will probably require substantially longer chains, much larger systems, and future computers. ACKNOWLEDGMENTS This work was supported by the Deutsche Forschungsgemeinschaft, and the Hochstleistungsrechenzentrum Jiilich within the projects "Polymer Dynamics" and "Thermodynamics of Disordered Polymer Systems." Stimulating discussions with C. Pierleoni, J.-P. Ryckaert, K. Binder, and G. S. Grest are gratefully acknowledged. lB. Diinweg and K. Kremer, Phys. Rev. Lett. 61, 2996 (1991). 2M. Bishop, M. H. Kalos, and H. L. Frisch, J. Chern. Phys. 70, 1299 (1979). 3D. C. Rapaport, J. Chern. Phys. 71, 3299 (1979). 4yU. Ya. Gotlib, N. K. Balabaev, A. A. Darinskii, and I. M. Neelov, Macromolecules 13, 602 (1980). 5W. Bruns and R. Bansal, J. Chern. Phys. 74, 2064 (1981). 6p. G. Khalatur, Yu. G. Papulov, and A. S. Pavlov, Mol. Phys. 58, 887 (1986). 7B. Smit, A. van der Put, C. J. Peters, J. de Swaan Arons, and J. P. J. Michels, J. Chern. Phys. 88, 3372 (1988). 8J. Luque, J. Santamaria, and J. J. Freire, J. Chern. Phys. 91, 584 (1989). 9C. Pierleoni and J.-P. Ryckaert, Phys. Rev. Lett. 61; 2992 (1991). IOC. Pierleoni and J.-P. Ryckaert, J. Chern. Phys. 96, 8539 (1992). llW. Smith and D. C. Rapaport, Mol. Sim. 9, 25 (1992). 12J. G. Kirkwood and J. Riseman, J. Chern. Phys. 16, 565 (1948). 13J. J. Erpenbeck and J. G. Kirkwood, J. Chern. Phys. 29, 909 (1958). 14B. H. Zimm, J. Chern. Phys. 24, 269 (1956). 15p. E. Rouse, J. Chern. Phys. 21, 1272 (1953).. 16M. Doi and S. F. Edwards, The Theory ofPolymer Dynamics (Clarendon, Oxford, 1986). 17D. L. Ermak and J. A. McCammon, J. Chern. Phys. 69, 1352 (1978). 18M. Fixman, Macromolecules 14, 1710 (1981). 19M. Fixman, J. Chern. Phys. 78,1594 (1983). 20M. Fixman, Macromolecules19, 1195 (1986). 21W. Zylka and H. C. Ottinger, J. Chern. Phys. 90, 474 (1989). 22A. R.ey, J. J. Freire, and 1. Garcia de la Torre, Macromolecules 24, 4666 (1991). 23K. Kremer and K. Binder, Comput. Phys. Rep. 7, 259 (1988). 24p. G. de Gennes, Scaling Concepts in Polymer Physics (Cornell University, Ithaca, 1979). 25B. H. Zimm, Macromolecules 13, 592 (1980). 26J. Batoulis and K. Kremer, Macromolecules 22, 4277 (1989). 27B. J. Berne and R. Pecora, Dynamic Light Scattering (Wiley, New York, 1976). 28y. ~Tsunashima, M. Hirata, N. Nemoto, and M. Kurata, Macromolecules 20, 1992 (1987). 29M. Bhatt, A. M. Jamieson, and R. G. Petschek, Macromolecules 22, 1374 (1989). 30B. Diinweg, J. Chern. Phys. 99, 6977 (1993), preceding paper. 31J. M. V. A. Koelman, Phys. Rev. Lett. 64,1915 (1990). 32p. J.Hoogerbrugie and J. M. V. A. Koelman, Europhys. Lett. 19, 155 (1992). 33K. Kremer, G. S. Grest, and I. Carmesin, Phys. Rev. Lett. 61, 566 (1988). 34K. Kremer and G. S. Grest, J. Chern. Phys. 92,5057 (1990). 35H. Risken, The Fokker-Planck Equation (Springer, Berlin, 1984). 36A. Z. Akcasu, M. Benmouna, and C. C. Han, Polymer 21, 866 (1980). 37M. Benmouna and A. Z. Akcasu, Macromolecules 13, 409 (1980). 38B. Diinweg and K. Kremer (to be published). 39M. O. Robbins, G. S. Grest, and K. Kremer, Phys. Rev. B 42, 5579 ( 1990). 4OB. R. A. Nijboer and F. W. de Wette, Physica 23,309 (1957). 41C. W. J. Beenakker, J. Chern. Phys. 85, 1581 (1986). 42G. S. Grest, B. Diinweg, and K. Kremer, Comput. Phys. Commun. 55, 269 (1989). _ J. Chern. Phys., Vol. 99, No.9, 1 November 1993 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions B. DOnweg and K. Kremer: Dynamics of a polymer chain 6997 43 J. D. Weeks, D. Chandler, and H. C. Andersen, J. Chern. Phys. 54, 5237 (1971). #T. Schneider and E. Stoll, Phys. Rev. B 17, 1302 (1978). 4SC. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations (Prentice-Hall, Englewood Cliffs, NJ, 1971). 046S. Toxvaerd, Mol. Phys. 56, 1017 (1985). 47J. P. Hansen and I. R. McDonald, Theory of Simple Liquids (Academic, New York, 1986). 48W. Paul, K. Binder, D. W. Heennann, and K. Kremer, J. Chern. Phys. 95,7726 (1991). 49C. Pierleoni (private communication). J. Chern. Phys., Vol. 99, No.9, 1 November 1993 Downloaded 17 Apr 2012 to 131.111.115.120. Redistribution subject to AIP license or copyright; see http://jcp.aip.org/about/rights_and_permissions