The glass transition of a hard sphere system is investigated within the
framework of the density functional theory (DFT). Molecular dynamics (MD)
simulations are performed to study dynamical behavior of the system on the one
hand and to provide the data to produce the density field for the DFT on the
other hand. Energy landscape analysis based on the DFT shows that there appears
a metastable (local) free energy minimum representing an amorphous state as the
density is increased. This state turns out to become stable, compared with the
uniform liquid, at some density, around which we also observe sharp slowing
down of the $alpha$ relaxation in MD simulations.Comment: 5 pages, 5 figure

A.D. Haymet

A.D.J. Haymet

A.R. Denton

B.D. Lubachevsky

B.J. Alder

C. Kaur

C.H. Bennett

D. Henderson

D.W. Oxtoby

H. Löwen

J. Jäckle

Kang Kim

M. Baus

M.D. Ediger

P. Tarazona

R. Yamamoto

T.V. Ramakrishnan

Toyonori Munakata

W. Kob

W.A. Curtin

X. Xia

Y. Rosenfeld

Y. Singh

English

arXiv

Crossref

Glass transition of hard sphere systems: Molecular dynamics and density functional theory

Kim, K

Munakata, T

Kyoto University Research Information Repository

Title Glass transition of hard sphere systems: Molecular dynamicsand density functional theoryAuthor(s)Kim, K; Munakata, TCitationPHYSICAL REVIEW E (2003), 68(2)Issue Date2003-08URL http://hdl.handle.net/2433/50290RightCopyright 2003 American Physical SocietyType Journal ArticleTextversionpublisherKyoto UniversityGlass transition of hard sphere systems: Molecular dynamics and density functional theoryKang Kim* and Toyonori Munakata†Department of Applied Mathematics and Physics, Graduate School of Informatics, Kyoto University, Kyoto 606-8501, Japan~Received 5 November 2002; published 12 August 2003!The glass transition of a hard sphere system is investigated within the framework of the density functionaltheory ~DFT!. Molecular dynamics ~MD! simulations are performed to study the dynamical behavior of thesystem on the one hand and to provide the data to produce the density field for the DFT on the other hand.Energy landscape analysis based on the DFT shows that there appears a metastable ~local! free energy mini-mum representing an amorphous state as the density is increased. This state turns out to become stable,compared with the uniform liquid, at some density around which we also observe a sharp slowing down of thea relaxation in the MD simulations.DOI: 10.1103/PhysRevE.68.021502 PACS number~s!: 64.70.Pf, 64.60.Cn, 61.20.JaUnderstanding the universal mechanism of the glass tran-sition is one of the major challenges in current condensedmatter physics. From a dynamical point of view, we wouldlike to know how a drastic slowing down near the transitionpoint ~temperature or density! and an eventual exceeding ofthe relaxation time over the experimental time scale are re-alized @1,2#. Energetically or statically, it is asked whether athermodynamic glassy state with an amorphous arrangementof particles has lower free energy than a liquid state of uni-form density. In order to answer these questions, many ef-forts have been devoted to real experiments, computer simu-lations, and theories in the last few decades.As one of the theories to study supercooled liquids andglasses, the density functional theory ~DFT! is recently gath-ering much attention @3#. The DFT is now a conventionalmethod to study the freezing @4,5# and other transitions. Theglass transition has been investigated also based on the DFTby some workers @6–10#. In the earlier works @6–10#, therandom close packing ~RCP! of hard spheres has been pro-duced by Bennett’s algorithm @11# and the free energy fromthe DFT with the input density field supplied by the RCPdata has been calculated. Singh et al. @6# showed that thehard sphere glassy state becomes more stable than a uniformliquid at the critical density ngs351.14, suggesting thatthere exists a kind of thermodynamic ~later called a randomfirst-order @12#! glass transition. Here, s is the hard spherediameter and n is the number density of the system. It isremarked here that since the RCP configurations were pro-duced by a kind of aggregation method, we cannot study thedynamic aspects of the glass transition found by the energet-ics based on DFT.The purpose of this paper is first to produce a supercooledand a glassy state for a one-component hard sphere system,relying on the uniform compression molecular dynamics~MD! method recently developed by Lubachevsky and Still-inger @13#, and then to study both dynamic and static prop-erties of the state. Especially from the particle configurationdata, we can discuss free energy within the DFT framework.This MD approach, in conjunction with the DFT, enables usto study both dynamic and static aspects of the transition incontrast with Bennett’s approach.Our system consists of N51372 identical hard sphereswith mass m and diameter s in a cubic box of volume V withperiodic boundary conditions. Throughout this paper, theunits of length and time are s and Ams2/kBT , respectively,where kB is Boltzmann’s constant and T is the temperature@14#. It should be noted that the temperature is fixed askBT51 in the course of the MD simulations.To begin with, we briefly explain our MD method to ob-tain glassy states of a one-component hard sphere system.Employing the standard Alder and Wainwright algorithm@15,16#, we first generate the equilibrium liquid state at den-sity n˜50.86. It is well known that the fluid system freezes atn˜ f.0.94. To avoid crystallization and to obtain amorphousglassy states, Lubachevsky and Stillinger introduced a com-pressing ~or quenching! procedure @13#, in which they actu-ally increased the diameter s with a constant rate of expan-sion during MD simulations. The dimensionless expandingrate G is defined asG5ds~ t !dt AmkBT, ~1!and G50.01 is chosen in our simulations. From the initialstate n˜50.86, we expanded each sphere with the rate G andcould obtain various high-density states n˜50.86, 0.94, 1.02,1.06, 1.10, 1.14, 1.18, and 1.21, without crystallization.Let us first study the static structure of the system. Forthis purpose, the radial distribution function g(r), which isdefined byg~r !51nN K (iÞ jNd~r1ri2rj!L , ~2!is calculated, where ri represents the positions of the ith par-ticle and ^& denotes the ensemble average over differentconfigurations. In Fig. 1, we plotted g(r) for densities n˜50.94, 1.06, 1.14, and 1.21. It is noted that there is no signof crystallization, which would be reflected in the sharp*Present address: Department of Physics, Kyoto University,Kyoto 606-8502, Japan. Electronic address: kin@scphys.kyoto-u.ac.jp†Electronic address: munakata@amp.i.kyoto-u.ac.jpPHYSICAL REVIEW E 68, 021502 ~2003!1063-651X/2003/68~2!/021502~5!/$20.00 ©2003 The American Physical Society68 021502-1peaks of g(r) at some characteristic lattice spacings. Instead,g(r) for higher density (n˜>1.14) shows splitting of the sec-ond peak, which is a familiar characteristic of a glassy stateof a simple liquid. Furthermore, we notice that the contactvalue g(r˜51) shows an anomalous increase as the densityincreases ~see the inset!, which corresponds to the fact thatthe pressure increases drastically with increasing n˜ . Inciden-tally, it should be remarked that the forms of g(r) agreequalitatively with those illustrated in Ref. @11#.We next consider the dynamic aspects of the hard sphereglasses by calculating the incoherent intermediate scatteringfunction Fs(q ,t), which is defined byFs~q ,t !5K 1N (j51Nexp@ iqDrj~ t;t0!#Lt0, ~3!where Drj(t;t0)5rj(t1t0)2rj(t0) is the displacement vec-tor of the j th particle in time t and ^& t0 represents anaverage over initial times t0. It is noted that Fs(q ,t) is one ofthe standard quantities in the studies of dynamic propertiesof supercooled liquids and glasses @17,18#. In Fig. 2, weplotted the decay profiles of Fs(q ,t) at a dimensionless wavenumber q˜52p for n˜50.94, 1.06, 1.14, and 1.21. We see inFig. 2 that the relaxation of Fs(q ,t) at n˜50.94 can be ex-pressed by a simple exponential function. Beyond the densityn˜51.06, however, Fs(q ,t) exhibits a two step, that is, fast band slow a , relaxation, which is often mentioned as a char-acteristic sign of the slow relaxation in glass forming liquids.At the highest density n˜51.21, the Fs(q ,t) does not showany decaying behavior @19#. One can define the structuralrelaxation time t˜ by Fs(q˜52p ,t˜ )5e21, and this t˜ is plot-ted as a function of n˜ in Fig. 3. We notice in Fig. 3 that therelaxation time t˜ shows a strong dependence on density,which can be expressed by the power law (n˜ g ,MD2n˜ )2gwith n˜ g ,MD.1.15 and g.1.31 ~solid line in Fig. 3!.We now consider energetics of the system based on theDFT and the configuration generated by the MD simulations.For a practical calculation of the DFT, we employ the Ra-makrishnan and Yussouff ~RY! free energy functional @4# be-cause of its simplicity and physical clarity. The RY func-tional is given byF@n#5Fid1Fint(2)1Funi , ~4!whereFid5kBTE n~r!lnFn~r!n Gdr, ~5!Fint(2)5212 kBTE E @n~r!2n#C~ ur2r8u!@n~r8!2n#drdr8.~6!Here, Fid and Funi represent the ideal gas contribution andthe excess free energy of the uniform liquid state n(r)5n ,respectively. Fint(2) represents the second-order term in the ex-pansion around the uniform liquid state, thus all terms higherFIG. 1. The radial distribution function g(r˜) obtained for n˜50.94 ~solid line!, 1.06 ~dashed line!, 1.14 ~short dashed line!, and1.21 ~dot-dashed line!. Inset: contact value of radial distributionfunction g(r˜51) as a function of n˜ . The units of r˜ and n˜ are s(r˜5r/s) and s23 (n˜5ns3), respectively.FIG. 2. Intermediate scattering function Fs(q˜ , t˜) at a wave num-ber q˜52p for n˜50.94 ~solid line!, 1.06 ~dashed line!, 1.14 ~shortdashed line!, and 1.21 ~dot-dashed line!. The units of q and t ares21 (q˜5qs) and Ams2/kBT ( t˜5tAkBT/ms2), respectively.FIG. 3. Structural relaxation time t˜ as a function of density n˜~closed circles!. Solid line represents power-law fit (n˜ g ,MD2n˜ )2gwith n˜ g ,MD51.15 and g51.31.K. KIM AND T. MUNAKATA PHYSICAL REVIEW E 68, 021502 ~2003!021502-2than second order are neglected. We note that C(r) is thedirect correlation function of the uniform liquid with densityn @20#.In order to evaluate the free energy of the system, weneed the trial density field n(r), for which we employ aconventional Gaussian superposition @6–10#; that is, the den-sity field n(r) is expressed by a sum of Gaussians with thecenters located at N sites $ri%, which are given by our MDsimulation,n~r!5(i51N S ap D3/2exp@2a~r2ri2!#[(i51Nz~r;ri!, ~7!where a21, a variational parameter for the calculation of thefree energy, is proportional to the mean square displacementof each particle. Small ~large! a represents the uniform liq-uid ~localized amorphous! state.When a is very large, Fid is asymptotically representedby @6#Fid~a!;NkBTF H 32 lnS ap D232J 2ln nG . ~8!For a small a region, we calculated the integral Eq. ~5! nu-merically. We confirmed that Fid approaches zero when a˜→0 and noticed that Fid coincides with the value of Eq. ~8!for a˜ *20.It is easy to see that the interaction term Fint can be di-vided into three parts as @6#Fint(2)~a!5212 NkBT H Fint ,1~a!1Fint ,2~a!2nE C~r !drJ ,~9!where Fint ,1(a) represents the self-interaction of a singleGaussian,Fint ,1~a!5E E z~r;0!C~ ur2r8u!z~r8;0!drdr8, ~10!and Fint ,2(a) represents the interaction between the two dis-tinct Gaussians,Fint ,2~a!5nE g~r1!dr1E E z~r;0!C~ ur2r8u!3z~r8;r1!drdr8. ~11!The pair distribution function g(r) in this equation is deter-mined from the MD simulation. With respect to the directcorrelation function C(r), we use Henderson-Grundke ex-pression for C(r), which is known to be reliable, thoughempirical, even for high-density hard sphere liquids @21#.The total free energy per particle relative to uniform state,D f ~a!5Fid~a!1Fint(2)~a!NkBT, ~12!is calculated as a function of the localization parameter a˜~see Fig. 4!. It is seen in Fig. 4 that the free energy localminimum at finite a , which represents an amorphous state,appears as the density is increased. As mentioned in Ref. @6#,the local minimum appears as the result of the competitionbetween the ideal gas ~simple increasing function of a) andthe interaction terms. Figures 4 show that two local minimaare located at a˜ .13 and 1600 for n˜51.14. Das and Kauralso observed two local minima of D f (a), which are calledthe weakly localized state for small a and the highly local-ized state for large a @9#. In addition, similar values for ahave been reported in Refs. @6,9#, which also use the RYfunctional, in relation to the local minimum of D f (a). Asstated in Ref. @9#, the qualitative adequacy is ambiguous forthe highly localized state with very high value of a since theRY form includes a perturbative expansion around the uni-form state. In fact, based on the MD data, we estimated a forhigh-density states by relating it to the plateau value of thetime-dependent mean square displacement of each particle,yielding a˜ .50 for n˜51.14, for instance. From this, it isseen that the RY form does not give the proper estimation ofthe degree of localization compared with the results obtainedin earlier works @8,10#.In Fig. 5, we plotted the free energy differences D f of theweakly and highly localized states as a function of density n˜ .From Fig. 5, we notice that the weakly localized state ap-pears for n˜>1.06 and the highly localized state appears forn˜>1.14. For higher densities n˜*1.15, it is seen that theFIG. 4. Total free energy per particle relative to uniform liquidD f (a˜ ) as a function of localization parameter a˜ for a˜ <100 ~a! anda˜ >100 ~b!. The unit of a˜ is s22 (a˜ 5as2).GLASS TRANSITION OF HARD SPHERE SYSTEMS: . . . PHYSICAL REVIEW E 68, 021502 ~2003!021502-3highly localized state is more stable than the weakly local-ized state. Moreover, it is found in Fig. 5 that both weaklyand highly localized states become more stable than the uni-form state at around n˜ g ,DFT.1.15, which is the liquid-glasstransition density from the energetics based on the DFT. Inpassing, we note that our ~random first order! glass transitiondensity n˜ g ,DFT51.15 is rather close to that n˜ g51.14, foundin Ref. @6#.Finally, we compare our results from the energetics abovewith dynamic information supplied by our MD simulations.We find in Fig. 2 that the intermediate scattering functionFs(q ,t) begins to exhibit the two-step relaxation at densityn˜.1.06, which corresponds precisely to the density wherethe free energy local minimum begins to appear in our DFT~see Figs. 4!. Turning to the relaxation time t˜ , we recall thatthe density dependence of t˜ could be described by the powerlaw (n˜ g ,MD2n˜ )2g with n˜ g ,MD.1.15. This density happenedto coincide with the density ng ,DFT , beyond which the local-ized state is more stable than the uniform liquid in thepresent DFT. From these results, we consider that the DFTbased energetics and dynamic behaviors related to slow dy-namics are well correlated with each other.In this paper, we reconsidered the DFT approach to theglass transition in the hard sphere system, which was firstundertaken by Singh et al. We obtained hard sphere glassesby MD simulations without recourse to the Bennett algo-rithm, and the information on particle configurations pro-duced by the MD simulations is used as input data when thefree energy is calculated based on the DFT. While only theuniform liquid state is stable at low density, the free energylocal minimum begins to appear at high density n˜.1.06,where our MD shows that two-step relaxation begins to ap-pear. This metastable glassy state becomes stable relative tothe uniform liquid at n˜ g ,MD51.15. Slow relaxation, as rep-resented by Fs(q ,t), turned out to be consistent with theenergetics based on the DFT.Before concluding this paper, we comment on recent de-velopments in the studies of the DFT. In recent years, theso-called weighted density approximation ~WDA! for thefree energy functional has been developed @22,23# and themodified version has also been introduced @24#. Moreover,the fundamental measure theory ~FMT! has been proposed@25# and is gathering considerable interest. It is well knownthat, for highly localized states, such methods are more ac-curate than the RY functional. This is because the formeremploys a nonperturbation approximation, whereas the latterrelies on a perturbation expansion around the uniform state.Several workers have already investigated the glass transi-tion by using the modified WDA method and found that themetastable localized state is located at a˜ .100 @8,10#, in ac-cordance also with our MD results. Although the DFT basedon the RY functional is still useful because of its physicalclarity, generality, and simplicity, we think that in view of therecent achievements, it is meaningful to employ the WDA orFMT method in order to obtain improved results.Furthermore, we expect that the present DFT approachwill be applied to more complex systems. As a model of aglass forming liquid, the binary Lennard-Jones @17# or soft-core system @18# has been investigated by large scale MDsimulations. Our approach can be readily applied to such asystem and would give new insights into the glass transitionfrom both thermodynamic and dynamic points of view.@1# J. Ja¨ckle, Rep. Prog. Phys. 49, 171 ~1986!.@2# M.D. Ediger, C.A. Angell, and S.R. Nagel, J. Phys. Chem. 100,13 200 ~1996!.@3# For reviews, see A.D.J. Haymet, Annu. Rev. Phys. Chem. 38,89 ~1987!; D.W. Oxtoby, in Liquid, Freezing, and the GlassTransition, edited by J.P. Hansen, D. Levesque, and J. Zinn-Justin ~Elsevier, New York, 1990!.@4# T.V. Ramakrishnan and M. Yussouff, Phys. Rev. B 19, 2775~1979!.@5# A.D. Haymet and D.W. Oxtoby, J. Chem. Phys. 74, 2559~1981!; D.W. Oxtoby and A.D. Haymet, ibid. 76, 6262 ~1982!.@6# Y. Singh, J.P. Stoessel, and P.G. Wolynes, Phys. Rev. Lett. 54,1059 ~1985!.@7# M. Baus and J.L. Colot, J. Phys. C 19, L135 ~1986!.@8# H. Lo¨wen, J. Phys.: Condens. Matter 2, 8477 ~1990!.@9# C. Kaur and S.P. Das, Phys. Rev. Lett. 86, 2062 ~2001!.@10# C. Kaur and S.P. Das, Phys. Rev. E 65, 026123 ~2002!.@11# C.H. Bennett, J. Appl. Phys. 43, 2727 ~1972!.@12# X. Xia and P.G. Wolynes, Proc. Natl. Acad. Sci. U.S.A. 97,2990 ~2000!.@13# B.D. Lubachevsky and F.H. Stillinger, J. Stat. Phys. 60, 561~1990!; B.D. Lubachevsky, F.H. Stillinger, and N. Pinson, ibid.60, 501 ~1991!.@14# For all physical quantities, the following reduced units areused: r˜5r/s , t˜5tAkBT/ms2, n˜5ns3(n[N/V), q˜5qs ,and a˜ 5as2.@15# B.J. Alder and T.E. Wainwright, J. Chem. Phys. 31, 459~1959!.@16# M.P. Allen and D.J. Tildesley, Computer Simulations of Liq-FIG. 5. Free energy differences D f of the weakly localized state~solid line! and the highly localized state ~dashed line! as a functionof density n˜ .K. KIM AND T. MUNAKATA PHYSICAL REVIEW E 68, 021502 ~2003!021502-4uids ~Clarendon Press, Oxford, 1987!.@17# W. Kob and H.C. Andersen, Phys. Rev. E 52, 4134 ~1995!.@18# R. Yamamoto and A. Onuki, Phys. Rev. E 58, 3515 ~1998!.@19# It is remarked here that due to frequent collisions under com-pressed situations, a long time simulation at high density isprohibitively difficult to perform.@20# J.P. Hansen and I.R. McDonald, Theory of Simple Liquids~Academic Press, London, 1986!.@21# D. Henderson and E.W. Grundke, J. Chem. Phys. 63, 601~1975!.@22# P. Tarazona, Phys. Rev. A 31, 2672 ~1985!.@23# W.A. Curtin and N.W. Ashcroft, Phys. Rev. A 32, 2909 ~1985!.@24# A.R. Denton and N.W. Ashcroft, Phys. Rev. A 39, 4701 ~1989!.@25# Y. Rosenfeld, Phys. Rev. Lett. 63, 980 ~1989!.GLASS TRANSITION OF HARD SPHERE SYSTEMS: . . . PHYSICAL REVIEW E 68, 021502 ~2003!021502-5

http://repository.kulib.kyoto-u.ac.jp/dspace/bitstream/2433/50290/1/e021502.pdf

Glass Transition of Hard Sphere Systems: Molecular Dynamics and Density
  Functional Theory

Glass Transition of Hard Sphere Systems: Molecular Dynamics and Density Functional Theory

Abstract

Similar works

Full text

Available Versions

Crossref

Kyoto University Research Information Repository