by molecular-dynamics simulation

F.Ould Kaddour 1,2and D.Levesque 11Laboratoire de Physique Th´e orique,Universit´e Paris Sud,Bˆa timent 210,91405Orsay,France 2Institut de Physique,Universit´e de Tlemcen,BP119Tlemcen 13000Alg´e rie (February 1,2008)Abstract By using the Kirkwood formula,the friction coefficient of a solvated Brown-ian particle is determined from the integration on time of the autocorrelation function of the force that the solvent exerts on this particle.Extensive molec-ular dynamics simulations show that above a definite size of the studied sys-tems the value of the integral defining the friction coefficient goes to a quasi constant value (a plateau)when the upper bound on time increases.The minimal value of the system size where the integral exhibits this asymptotic behavior,rises with the Brownian particle size.From the plateau,a reliable estimate of the friction coefficient is obtained.PACS numbers:83.10.Mj,83.10.Rs In solutions,it is supposed that the large particles such as micelles or colloids which coexist with the atoms,ions or small molecules of the solvent behave as Brownian parti-

cles.At low concentrations of the large particles,by using multi-scale analysis [1–4],this hypothesis has been justified in the limit where the ratio between the mass of the solvated particles and that of the solvent molecules goes to ∞.It has then been established that the diffusion coefficient of Brownian particles can be computed in term of the friction coefficient characterizing the force exerted on them by the solvent.When the Brownian particles have

However in many suspensions,such as ionic solutions,the values of ratios of masses and sizes between the solvated particles and solvent molecules are only of the order of10and the possibility that the solvated particles can be considered as Brownian particles becomes questionable.Several theoretical works[6,7]and works based on numerical simulations [8–11]have been devoted to this question.The main problem addressed in these works was that of the determination of the lower bounds of the size and mass ratios above which,to a good approximation,the motion of the solvated particles is Brownian.The criterion chosen to locate these bounds was that the diffusion coefficient of the solvated particles obeys to the relation between the diffusion coefficient D and friction coefficientξstrictly valid only for brownian particles D=k B T/ξ(k B is Boltzmann’s constant,T the temperature of the solvent).The main concern,when the Stokes estimate ofξis used,is the choice of the hydrodynamic boundary condition between solvated particles and solvent which is well defined only when the solvated particle has a macroscopic size.This last shortcoming can be overcome by computingξfrom its exact expression for a Brownian particle derived by Kirkwood[12]and later,more rigourously,from multi-scale analysis[2,13].Obviously this method seems the correct way to proceed in order to check the brownian behavior of a solvated particle.However,as it was discussed in the literature[8,10,12,14–16],this method is not easy to use in simulations due to importantfinite size effects.This work is devoted to discuss this problem and to establish in what conditions,in a numerical simulation,the friction coefficient of a solvated particle of large mass and size can be credibly determined.

The friction coefficientξis given in terms of the integration on time t of the equilib-rium autocorrelation functionof the instantaneous microscopic force F(t) experienced by the Brownian particle:ξ=

1

dt ≡−˙P(t)=−

N

i=1

m

d v i(t)

3k B T lim

t→∞

t

N dτ

=1

dτ

1

3k B T

+

<˙P(0).P(0)>NIf in simulations N is large enough,it can be expected,following the remark made by Kirkwood in[12],that,in the range of values of t whereξN(t)reaches its asymptotic form≃c g(at/N),t is such as t<As it was proposed,for instance,in[6,10]it is possible to give a specific analytic form to g(at/N)by supposing that,following the Onsager’s principle,the regression of thefluctua-tions of F(t)at large t is governed by the laws of the macroscopic hydrodynamics.According to these laws,the force exerted by the solvent on the Brownian particle is proportional to the momentum P(t)carried by the solvent,i.e.

F(t)=

ξo

Nm

t)(6)

and then

ξN(t)=−<˙P(t).P(0)>N

Nm

t).(7)

At large N and t such as t<<ξo/Nm,an expansion of the exponential yields to a linear expression,similar to that ofξN(t)given above,allowing to determine the friction coefficient as

ξN(t)=ξo(1−

ξolarge compared of those of the solvent molecules.Hence the mass of this particle satisfies to the condition required so that the relation between D andξ,given above,applies.When M is large,it is possible to consider that,in the time scale accessible in a simulation,the particle is immobile.

The molecules and thefixed particle interact through a Lennard-Jones(LJ)potential modified with a cubic spline,as describe in a previous paper[11].This potential has the form v ij(r)=ǫij f(r/σij)where the indices i,j=1or2refer to the solvent andfixed particle,respectively.The parametersσij are such asσ12=(σ11+σ22)/2,andǫij are equal ǫ12=ǫ22=ǫ11.The unit of time is chosen equal toτ0=fit ofξN(t),for t/τ0between12and20,to an exponential form when N is≃1000and a linear form when N is≃20000.

However,the possibility that the values ofξ,determined from thefits made at different values of N,coincide within the statistical uncertainties,supposes thatfinite size effects,in particular those associated to the use of the periodic boundary conditions,do not affect the long time behavior ofξN(t).From thefit of the multiplicative constant of the exponential (cf.Eq.(7))the estimates ofξare:219.4(N=8),144.7(N=1500)and122.9 (N=5324)and from that of the constant term in the linear form(cf.Eq.8)they are: 113.3(N=12000)and110.7(N=32000).Clearly,for the two large systems these values of ξagree within the statistical error equal to10-15%.But the factor of2between the values found at N=8and N=32000must be attributed tofinite size effects.The side lenghts L of the simulation cells being for these two values of N equal to∼10and∼35,due to the periodic boundary conditions thefixed particle is distant from these nearest replica by the same lengths.The sound velocity c s of the solvent,for the considered thermodynamic state,being in reduced units∼6,gives typical times L/c s of1.5and6beyond of which thefinite size effects resulting from the mutual influence between thefixed particle and its replicas can affect the correlation function.The magnitude of these effects on the value ofξis difficult to assess.It has been quantitatively discussed only for the bulk values of the transport coefficients of diffusion or viscosity in densefluids[18],corrections of about ∼10%have been found for systems of N∼1000molecules and they seem much larger on ξ.Other estimates ofξare obtained from thefit of the coefficient of t in the exponential (Eq.(7))or the linear approximation(Eq.(8)).For the small values N,we obtained205.5 (N=8)and143.3(N=1500),for the large values of N,98.4and91.1.These results seem to confirm that the asymptotic behavior ofξN(t)is well described by Eq.(7)taking into account thefinite size effects.

FIGURES

t/τ0050

100

150

200

ξΝ(t )

FIG.1.ξN (t )for the fixed particle of size σ22=4.0and increasing values of N .Insert :asymptotic behavior with error bars and its linear fit at N =12000and 32000

t/τ00100

200

300400

500

ξΝ(t )FIG.2.ξN (t )for the fixed particle of size σ22=7.0and increasing values of N .Insert :asymptotic behavior wih error bars and its linear fit for N =32000and 55296

We consider now the case of the fixed particle with a size σ22=7.0.The ξN (t )functions are plotted in Fig.2.for N equal to 5324,12000,32000and 55296.The behavior of ξN (t ),when N increases,is similar to that obtained with the particle of size σ22=4.0,but it is only for N ≥32000that ξN (t )presents a slow linear decrease for t/τ0>10.0.For the particle of size σ22=4.0,we remark that this latter behavior is reached for systems smaller by a factor ∼2and,then,the size of the fixed particle has an important influence on the value of the system size where ξN (t )exhibits a slow linear decrease at large time.This remark is confirmed by the results of the fits of ξN (t )by an exponential at N =5324and 12000or a linear function at N =32000and 55296.The values of ξare 217.0,244.2,209.3and 207.5respectively.We notice that,the linear behavior of ξN (t ),at large t ,corresponds obvioulsy to the fact that decreases very slowly in the same domain of t .This point is illustrated in Fig.3where this behavior is clearly seen for N =12000,32000and 55296within statistical errors.

0510

15t/τ0−700−600

−500

−400

−300

−200

−100

N 051015t/τ0051015

t/τ

0FIG.3.Autocorrelation function N for the particle of size σ22=7.0at N =12000,32000and 56296.Solid lines :estimate of its almost constant value from its av-erage for t/τ0>12.0and <20.0:∼-12.6,-6.6and -6.1.

From the results of the present simulations it seems needed to adopt a critical point of

As mentionned already,the asymptotic form ofξN(t)at large times has been discussed in many works in the literature,for instance in[14],[15]and[16],in particular the question of the occurence of a domain of time where the friction coefficientξN(t)should exhibit a plateau. It has been proposed in[8]to bypass the search of such a plateau inξN(t)by computing ξfrom the integration offrom t=0to the value of t=t1where,for the first times when t increases,becomes0.In our simulations t1/τ0≃0.5, clearly from the comparison between Fig.2and Fig.3ofξN(t)andN,such a method to estimateξseems problematic since the value ofξN(t1)for instance at N=55296 does not agree with the valueξobtained from the analysis of the asymptotic behavior of ξN(t).

Since the present simulations show that system sizes of N≃20000are needed to correctly estimatedξ,it is expected that similar system sizes are needed to compute D in order to avoidfinite size effects.For instance in[11]for N=5324,for a brownian particle with σ22=4.0and M=60in a LJ solvent at a thermodynamic state identical to that considered here,it was found D=0.077.By using D=k B T/ξthis value of D agrees well with that of ξ≃120obtained in this work at N=5324,but it differs by10%from that,computed atN=32000,ξ≃110.A new determination of D at N=32000gives D=0.095,which now agrees with our estimate ofξat the same value of N from the slow linear decay ofξN(t)for t/τ0>12.

In conclusion from our simulations realized for increasing system sizes from N∼1000to ∼56000,we have shown that the qualitative behavior ofand,consequently, that ofξN(t)is in excellent agreement with that guessed by Kirwood[12]and in subsequent works.[15,16].For large times,the representation ofξN(t)by Eq.(7),derived from a simple argument based on Onsager principle,seems an adequate model of this asymptotic behavior. From a quantitative point of view,in spite of simulation runs totalizing4to8millions of time steps the statistical uncertainties on the friction coefficient stay of the odrer of15%,a reduction of this uncertainty by a order of magnitude seems beyond what it is possible to make by using present standard computers.We stress that the aim of this work was the investigation of the variation ofξN(t)with N.In the thermodynamic limit,the asymptotic behavior ofξN(t)should be an algebraic decay at very large time.It was not considered here,in particular,because its amplitude is much smaller than the present uncertainties on ξN(t)in this domain of time[11,18].REFERENCES

[1]E.Frieman,Jour.Math.Phys.4,410(1963).

[2]R.I.Cukier and J.M.Deutch,Phys.Rev.177240(19).

[3]L.Bocquet,J.Piasecki and J-P Hansen,J.Stat.Phys.76505(1994).

[4]J.Piasecki,L.Bocquet and J-P Hansen,Physica A218125(1995).

[5]L.D.Landau and E.M.Lifshitz,Fluid Mechanics(Pergamon Press,London,1963).

[6]L.Bocquet,J-P Hansen and J.Piasecki,J.Stat.Phys.76527(1994).

[7]S.Bhattacharyya and B.Bagchi,J.Chem.Phys.106,1757(1997).

[8]A.N.Largar’kov and V.H.Sergeev,Sov.Phys.Usp.27,566(1978).

[9]J.J.Brey and J.G´o mez Ord´o˜n ez,J.Chem.Phys.76,3260(1982).

[10]P.Espa˜n ol and I.Z´u˜n iga,J.Chem.Phys.98,574(1993).

[11]F.Ould-Kaddour and D.Levesque,Phys.Rev.E6,61,(2000).

[12]J.Kirkwood,J.Chem.Phys.14,180(1946).

[13]P.R´e sibois and M.De Leener,Classical kinetic theory offluids(John Wiley and Sons,

Inc,1977).

[14]A.Suddaby and P.Gray,Proc.Phys.Soc.A75,,109(1960).

[15]E.Helfand,Phys.Fluids4681(1961).

[16]H.Mori,Prog.Theo.Phys.33,423(1965).

[17]M.Allen and D.Tildesley,Computer Simulation of Liquids(Oxford University Press,

Oxford,England,1987).

[18]J.Erpenbeck,Phys.Rev.A38,6255(1988).下载本文