A new Monte Carlo algorithm is introduced for the simulation of supercooled
liquids and glass formers, and tested in two model glasses. The algorithm is
shown to thermalize well below the Mode Coupling temperature and to outperform
other optimized Monte Carlo methods. Using the algorithm, we obtain finite size
effects in the specific heat. This effect points to the existence of a large
correlation length measurable in equal time correlation functions.Comment: Proceedings of "X International workshop on Disordered Systems" held
  in Molveno (Italy), March 200

Fernandez, L. A.

Martin-Mayor, V.

Verrocchio, P.

English

arXiv

© 2007 Taylor & Francis. International Workshop on Disordered Systems Moveno-Andalo (10. 2006. Trento, Italy). We thank G. Biroli, G. Parisi and R. de Nalda for discussions. We were partly supported by BSCH-UCM, CAM-UCM and by MEC (Spain), through contracts BFM2003-08532, FIS2004-05073 and FPA2004-02602. Simulations were carried out at BIFI.A new Monte Carlo algorithm is introduced for the simulation of supercooled liquids and glass formers, and tested in two model glasses. The algorithm thermalizes well below the Mode Coupling temperature and outperforms other optimized Monte Carlo methods.BSCH-UCM (Spain)CAM-UCM  (Spain)MEC (Spain)Depto. de Física TeóricaFac. de Ciencias FísicasTRUEpu

Fernández Pérez, Luis Antonio

Martín Mayor, Víctor

Docta Complutense

For Peer Review Only       Optimized Monte Carlo Method for glasses   Journal: Philosophical Magazine & Philosophical Magazine Letters Manuscript ID: TPHM-06-May-0145.R2 Journal Selection: Philosophical Magazine Date Submitted by the Author: 10-Jul-2006 Complete List of Authors: Fernandez, Luis; Universidad Complutense de Madrid, Departamento de Fisica Teorica Martin-Mayor, Victor; Universidad Complutense de Madrid, Departamento de Fisica Teorica Verrocchio, Paolo; University of Trento, Department of Physics Keywords: glass transition, Monte-Carlo Keywords (user supplied):     Note: The following files were submitted by the author for peer review, but cannot be converted to PDF. You must view these files (e.g. movies) online. molveno_revised2.tex molveno_revised2.bbl 3hist-global-local.epsi fase_desordenada.epsi swap.epsi swap_vs_micro.epsi    http://mc.manuscriptcentral.com/pm-pmlPhilosophical Magazine & Philosophical Magazine LettersFor Peer Review OnlyOptimized Monte Carlo Method for glassesL. A. Fernández1,3,∗, V. Mart́ın-Mayor1,3,† and P. Verrocchio2,3,4‡1 Dep. de F́ısica Teórica I, U. Complutense, 28040 Madrid, Spain.2 Dip. di Fisica, U. di Trento, 38050 Povo, Trento, Italy.3Instituto de Biocomputación y F́ısica de Sistemas Complejos (BIFI). Zaragoza, Spain.4Research center Soft INFM-CNR c/o U. di Trento 38050 Povo, Trento, Italy.July 10, 2006AbstractA new Monte Carlo algorithm is introduced for the simulation of supercooled liquids and glassformers, and tested in two model glasses. The algorithm thermalizes well below the Mode Couplingtemperature and outperforms other optimized Monte Carlo methods.∗laf@lattice.fis.ucm.es, †victor@lattice.fis.ucm.es, ‡paolo.verrocchio@unitn.it1Page 1 of 6http://mc.manuscriptcentral.com/pm-pmlPhilosophical Magazine & Philosophical Magazine Letters123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960For Peer Review Only1. IntroductionThe lack of structural or thermodynamic changes at the glass transition [1] is a major problem forits investigation. The only standard feature, as compared with second-order phase transitions [2], isthe dramatic dynamical slowing down at the critical temperature. A fairly standard mechanism forslow dynamics in a homogeneous system at finite temperature is the divergence of a correlation length(critical slowing down [2]). Slowness arises from the need of configurational changes to propagate overincreasingly large regions (the critical origin of the Mode Coupling singularity[3] has been recentlyrecognized[4]). It has been recently proposed [5] to study this growing length-scale in glass-formersthrough the finite size behaviour of small systems [2]. Note that experiments in films and nanopores[6, 7]show that the glass transition changes in samples with one or more dimensions of nanometric scale. Inparticular, the specific heat is most sensitive to the size of the confining pore when temperatures areclose to the glass transition[8].Numerical simulations are an important tool for the study of the glass transition. Their worsedrawback is the shortness of the times that may be simulated in today computers (roughly speaking,microseconds). As a consequence, the computer model goes out of equilibrium by the Mode-Couplingtemperature, Tmc, rather than the actual glass temperature, Tc. To approach Tc, one may resort tooptimized Monte Carlo (MC) methods [9, 10, 5], namely methods implementing unphysical dynamicalrules that strongly reduce the equilibration times. When thermalization is achieved, optimized MCallow one to study equilibrium mean values and their temperature (or pressure) derivatives, althoughthe purely dynamic features of these methods are interesting in their own right [5].In this conference note, we give the first full description of the local swap algorithm [5]. We com-pare the performance of the local swap dynamics with the standard MC and with the microcanonicalalgorithm [10]. We conclude that local swap yield equilibrium data at temperatures where the micro-canonical algorithm no longer thermalizes. Finally, we address the fishy issue of estimating the specificheat in a metastable liquid state. We give here details about the strategy followed in [5], where tinybut clearly measurable finite-size effects were observed in the specific heat.2. Models and observablesWe consider two similar models of fragile glass formers, namely binary mixtures of soft spheres. Thefirst model, extensively studied in Ref. [5], is a 50% mixture of particles interacting through the pairpotential Vαβ(r) = ε[(σα + σβ)/r]12 + Cαβ , where α, β = A,B, with a cutoff at rc =√3σ0. The choiceσA = 1.2σB hampers crystallization, as compared with the σA = σB model. We impose (2σA)3+2(σA+σB)3 + (2σB)3 = 4σ30 where σ0 is the unit length. Constants Cαβ are chosen to ensure continuity atrc. We work at constant volume, with particle density fixed to σ−30 and temperatures in the range[0.897Tmc, 10.792Tmc]. We use periodic boundary conditions in a cubic box in systems with N =512, 1024, 2048 and 4096 particles. For argon parameters, σ0 = 3.4Å, ε/kB = 120K and Tmc = 26.4K.In the second model [11, 10] the choice σB = 1.4σA is made. Naming x = r/(σα + σβ), the pairpotential in units ε = 1 is V (x) = x−12+x−13/1212/13 if x < 121/13 and zero otherwise (thus, the cut-offdistance depends on the type of interacting species). We study density 1.08σ−30 , as in Refs. [11, 10].Since the potential energy per particle, e shows the slowest excitations [5, 11], we shall focus hereonly in this observable, the internal energy being 32kBT + 〈e〉 (for other quantities, see Ref. [5]). Theconstant-volume specific heat, Cv , is:e =12N∑j,k 6=jV (~rk − ~rj) , Cv =32+NT 2[〈e2〉 − 〈e〉2](units ε = kB = 1). (1)2Page 2 of 6http://mc.manuscriptcentral.com/pm-pmlPhilosophical Magazine & Philosophical Magazine Letters123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960For Peer Review OnlyFigure 1: The local swap move. The A particles are depicted by full symbols, while B particles areopen symbols. The (say) B particle picked randomly (signaled by a central mark) may be swappedwith any one of the three A particles inside the sphere shown in the left part of the plot. After theswap (right), only two B particles could be exchanged with the picked B particle.3. The local swap Monte Carlo algorithmThe Grigera-Parisi swap algorithm[9] consists in picking randomly a pair of particles of distinct type, Aand B, try to exchange their positions, and accepting this move with probability min{1, e−∆E/T} (∆Eis the total potential energy change produced by the swap). If combined with standard MC, it is veryeffective in reducing the equilibration time. The caveat is that the acceptance of the swap move is verysmall (see below, and Fig. 4—right). Actually, the probability of accepting a swap move decreases withthe distance of the swapped particles, being almost asymptotic already at rc. For the nearest pickedparticles, at T = Tmc, the probability of accepting the swap is about 10 times larger than for distantparticles. In large systems, one seldom picks neighboring particles for the swap. To cure this problem,we have proposed the local swap.In the local swap, the elementary MC step is either (with probability p) a single-particle displacementattempt or (with probability 1 − p) an attempt to swap particles. Therefore, for p = 1 the algorithmreduces to standard MC. From here on we call local swap to the algorithm with p = 0.5. The timeunit t0 is N/p elementary steps. Our swap consists in picking a particle at random (the picked particle)and trying to interchange its position with that of a particle of opposite type (the swapped particle),chosen at random among those at distance smaller than 0.6rc for the σA/σB = 1.2 model and 0.73rcfor the σA/σB = 1.4 one (see Fig. 1). Note that this distance may not decrease arbitrarily, since noswappable particles would be eventually found. The optimal choice for the swap distance is the onethat maximize the number of actually performed swap moves in a given simulation time. Yet, there isan intrinsic lack of symmetry. Indeed, let Nold be the number of swappable particles around the pickedparticle in its original position, and Nnew the number of swappable particles for the picked particle in itsfinal position. The probability of choosing the swapped particle is 1/Nold in the original configuration,while it would be 1/Nnew in the final configuration. Detailed balance (see e.g. [2]) holds only if thisasymmetry is incorporated in the probability of accepting the swap:paccept swap = min{1,NoldNnewe−∆E/T } . (2)The performance of the local swap below Tmc is far superior to the standard MC (Fig. 2) andto the microcanonical method of Ref. [10] (Fig. 3). Furthermore, for both the σA/σB = 1.2 andthe σA/σB = 1.4 models, local swap finds a crystallization phase transition. A previous simulationof the σA/σB = 1.2 model [12], using the Grigera-Parisi global swap, found crystals with up to 60particles. The local swap, (Fig. 4—left), does find crystals even for systems with 1024 particles. For theσA/σB = 1.4 model, we find disordered bcc crystals, not reported in previous studies [11, 10]. Actually,3Page 3 of 6http://mc.manuscriptcentral.com/pm-pmlPhilosophical Magazine & Philosophical Magazine Letters123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960For Peer Review Only 2.01 2.02 2.03 2.045×106 1×107 2×107et/t05×106 1×107 2×107t/t0Figure 2: Comparison of the performance of the local swap (right) and standard Monte Carlo (left),as shown by the Monte Carlo history of the potential energy per particle (data points are the averageof 104t0 successive steps), for the σA/σB = 1.2 model at T = 0.897Tmc and 1024 particles. The twostandard Monte Carlo runs had as starting points thermalized configurations obtained with local swap.In our time window, the two standard simulations do not explore the same energy range. Instead thesingle local swap simulation explores the full energy range.with 54 particles, we found a perfect bcc crystal (the A and B particles formed two interleaving singlecubic lattices).The acceptances of the Grigera-Parisi global swap and the local swap algorithms in the fluid phaseare compared for the σA/σB = 1.4 model in Fig. 4—right, for aN = 32 and a N = 256 system. We havechecked that, for N = 256 particles, the acceptance for the global swap is already in the thermodynamiclimit (the acceptance of the local swap is N independent). For this model, the acceptance probabilitydecreases by three orders of magnitude with decreasing temperature, in the plotted range. In thethermodynamic limit, local swap outperforms the global swap by a factor 4 at 2.13 Tc, that improves toa gain factor of 8 at 0.9 Tc. For the model σA/σB = 1.2, the local swap acceptance varies linearly in T−4from 0.74% at 0.9Tmc up to 6% at 2Tmc. For the global swap, the acceptance in the thermodynamiclimit for the σA/σB = 1.2 model range from 0.53% at 0.9 Tmc up to 4.7% at 2Tmc.4. Studying the metastable phaseThe liquid is metastable in our soft-sphere models. Thus, the question arises of how to study the ther-modynamic properties of a supercooled liquid. We give here details about the strategy we followed inRef. [5]. The underlying assumption is that the equilibration time (specifically, the exponential autocor-relation time for the energy [5]) for the metastable liquid phase is much smaller than the crystallizationtime. Our strategy has been to run several MC runs (up to 400 at the lowest temperatures [5]). In eachrun, the equilibrated metastable liquid is neatly separated from the crystallizing system (see Fig. 4-left).In the analysis we only consider histories whose metastable liquid part was self-consistently found tobe longer than 100 exponential energy autocorrelation times (to avoid bias, we also discarded some 20autocorrelation times before crystallization). Note that the central MC history in Fig. 4-left does notmeet this criterion (at this low temperature the autocorrelation time was 105t0).4Page 4 of 6http://mc.manuscriptcentral.com/pm-pmlPhilosophical Magazine & Philosophical Magazine Letters123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960For Peer Review Only 2.5 3 3.5 4 4.5 0  0.5  1  1.5  2CvT/Tmcmicrocanonicallocal swapFigure 3: Specific heat of the σA/σB = 1.4 model with 128 particles, as a function of temperature, asobtained with local swap (109t0 steps) and with the microcanonical method [10] (error estimates werenot provided in Ref. [10]). Local swap produces a significantly larger estimate of the specific heat atour lowest simulated temperatures, signaling better sampling of configuration space and indicating thatthe microcanonical method is unable to thermalize at such low temperatures.1.952.002.05 1×107  2×107  3×107et/t01.952.002.05e1.952.002.05e10-610-510-410-31.0 1.5 2.0T/TmcAcceptanceFigure 4: (Left) Examples from the 100 generated (local swap) Monte Carlo histories for the potentialenergy density of the σA/σB = 1.2 model with 1024 particles at T = 0.897Tmc. One may clearlydistinguish the metastable liquid from the crystallizing system, that starts as a sharp energy drop in allthree runs. (Right) Swap probability acceptance for the here explained local algorithm (squares) andthe global Parisi Grigera swap (circles), for the σA/σB = 1.4 model in the fluid phase, as a functionof temperature. We show results for N = 32 particles (filled symbols) and N = 256 particles (opensymbols). The latter coincides within errors with those for N = 1024 particles.5Page 5 of 6http://mc.manuscriptcentral.com/pm-pmlPhilosophical Magazine & Philosophical Magazine Letters123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960For Peer Review Only5. ConclusionsWe give the first full description of the local swap algorithm [5], that is compared with the Grigera-Parisi global swap [9]. The swap acceptance for the local swap algorithm is neatly superior. Wecompare as well the performance of the local swap dynamics with the standard MC and with themicrocanonical algorithm [10], concluding that local swap yield equilibrium data at temperatures wherethe microcanonical algorithm no longer thermalizes. Furthermore, local swaps finds crystal states notreported in previous work for the σA/σB = 1.4 model [10, 11]. We have given details about how toestimate the specific heat in a metastable liquid state [5].AcknowledgmentsWe thank G. Biroli, G. Parisi and R. de Nalda for discussions. We were partly supported by BSCH—UCM, CAM—UCM and by MEC (Spain), through contracts BFM2003-08532, FIS2004-05073 andFPA2004-02602. Simulations were carried out at BIFI.References[1] Debenedetti, P. G. (1997) Metastable Liquids, Princeton University Press, .[2] Amit, D. and Mart́ın-Mayor, V. (2005) Field Theory, the Renormalization Group and CriticalPhenomena, World Scientific Singapore, in press.[3] Götze, W. and Sjögren, L. (1992) Rep. Prog. Phys. 55, 241–336.[4] Biroli, G. and Bouchaud, J.-P. (2004) Europhys. Lett. 67, 21–27.[5] Fernandez, L. A., Martin-Mayor, V., and Verrocchio, P. (2006) Phys. Rev. E 73, 030501.[6] Alcoutlabi, M. and McKenna, G. B. (2005) J. Phys.: Condens. Matter 17, 461–524.[7] Alba-Simionesco, C., Coasne, B., Dosseh, G., Dudziak, G., Gubbins, K., Radhakrishnan, R., andSliwinska-Bartkowiak, M. (2006) J. Phys.: Condens. Matter 18, 15–68.[8] Morineau, D., Xia, Y., and Alba-Simionesco, C. (2002) J. Chem. Phys. 117, 8966–8972.[9] Grigera, T. S. and Parisi, G. (2001) Phys. Rev. E 63, 045102–1–045102–4.[10] Yan, Q., Jain, T., and dePablo, J. (2004) Phys. Rev. Lett. 92, 235701.[11] Yu, C. C. and Carruzzo, H. M. (2004) Phys. Rev. E 69, 051201–1–051201–10.[12] Brumer, Y. and Reichman , D. R. (2004) J. Phys. Chem. B 108. 6832.6Page 6 of 6http://mc.manuscriptcentral.com/pm-pmlPhilosophical Magazine & Philosophical Magazine Letters123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960

The Philosophical Magazine A Journal of Theoretical Experimental and Applied Physics

Optimized Monte Carlo method for glasses

L. A. Fernández

V. Martín-Mayor

P. Verrocchio

Crossref

International audienceA new Monte Carlo algorithm is introduced for the simulation of supercooled liquids and glass formers, and tested in two model glasses. The algorithm thermalizes well below the Mode Coupling temperature and outperforms other optimized Monte Carlo methods

Optimized Monte Carlo Method for glasses

Abstract

Similar works

Full text

Available Versions

Docta Complutense

Crossref

HAL: Hyper Article en Ligne