Population Annealing Simulations of a Binary Hard Sphere Mixture
PPopulation Annealing Simulations of a Binary Hard Sphere Mixture
Jared Callaham ∗ and Jonathan Machta
1, 2, † Department of Physics, University of Massachusetts, Amherst, Massachusetts 01003 USA Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, New Mexico 87501 USA
Population annealing is a sequential Monte Carlo scheme well-suited to simulating equilibriumstates of systems with rough free energy landscapes. Here we use population annealing to study abinary mixture of hard spheres. Population annealing is a parallel version of simulated annealingwith an extra resampling step that ensures that a population of replicas of the system representsthe equilibrium ensemble at every packing fraction in an annealing schedule. The algorithm and itsequilibration properties are described and results are presented for a glass-forming fluid composedof a 50/50 mixture of hard spheres with diameter ratio of 1.4:1. For this system, we obtain preciseresults for the equation of state in the glassy regime up to packing fractions ϕ ≈ .
60 and studydeviations from the BMCSL equation of state. For higher packing fractions, the algorithm falls outof equilibrium and a free volume fit predicts jamming at packing fraction ϕ ≈ . PACS numbers:
I. INTRODUCTION
One of the grand challenges of computational statisti-cal physics is to understand the nature of glassy systems.Profound open questions remain concerning both config-urational glasses and spin glasses. The signature prop-erty of glassy systems is an extreme slowing of dynamicsin both laboratory and computational experiments. Theequilibrium properties of glassy systems are thus verydifficult to study. The situation is somewhat better forspin glasses where it is established that there is ther-modynamic glass transition and a well-understand meanfield theory of the low temperature glass phase, thoughcontroversies surround the nature of the low temperaturephase in finite dimensions. For configurational glasses itis not even known whether there is an underlying ther-modynamic transition or whether the glass transition isentirely a kinetic phenomenon.Much of the progress for finite-dimensional spin glasseshas been made possible due to algorithmic advances,particularly the introduction of replica exchange MonteCarlo, also known as parallel tempering [1–5]. The situ-ation for simulating equilibrium glassy fluids is less well-developed and fundamentally a more difficult problem.Parallel tempering has also been extensively applied tofluid systems [6] in the glassy regime [7, 8]. Other re-cent algorithmic advances such as event chain MonteCarlo [9, 10] and particle-swap Monte Carlo [11] havealso shown promises.In this paper we introduce population annealing MonteCarlo as a method for studying fluid systems and showthat it is an effective tool for studying equilibrium prop-erties up to high densities in the glassy regime. Popula- ∗ Electronic address: [email protected] † Electronic address: [email protected] tion annealing was first developed for spin glasses [12–14]and has been shown to be an efficient method for large-scale studies of equilibrium states [15, 16] and groundstates [17] of spin glasses. Here we use population an-nealing Monte Carlo to simulate a binary mixture of hardspheres.In particular, we study a glass-forming 50/50 binarymixture of hard spheres with a diameter ratio of 1.4:1,which has been the subject of previous computationalstudies, e.g. [8, 18, 19]. Using population annealing weobtain precise estimates of the equation of state up topacking fractions in the glassy regime. We compareour results to previous simulations and the well-knownBoublik-Mansoori-Carnahan-Starling-Leland equation ofstate [20, 21].The paper is organized as follows. In Sec. II we intro-duce population annealing and describe its properties.In Sec. III we describe the hard sphere fluid model, thesimulation methods and the observables that we study.Results for the equation of state, the approach to jam-ming and the performance of the algorithm algorithm arepresented in Sec. IV. The paper concludes with a discus-sion.
II. POPULATION ANNEALING MONTECARLOA. Overview
Population annealing (PA) is closely related to simu-lated annealing. In both algorithms, a system is takenthrough an annealing schedule in one or more thermody-namic control parameters, e.g. temperature or density,from a region where equilibration is easy, e.g. high tem-perature or low density, to a region where equilibrationis difficult. At each step in the annealing schedule, anequilibrating procedure, e.g. the Metropolis algorithm, a r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l is applied to the system at the current value of the con-trol parameter(s). The objective of simulated annealingis to find ground states or low-lying states and simu-lated annealing does not sample from the equilibriumdistribution at each step along the annealing schedule.Population annealing, by contrast, does sample from theequilibrium distribution at each step along the anneal-ing schedule. In PA a large population of replicas of thesystem is annealed and, at each step in the annealingschedule, the population is resampled to maintain equi-librium as described in the next section. B. Population Annealing for Hard Spheres
In this section we describe population annealing forhard sphere fluid systems using packing fraction ϕ as thecontrol parameter. A population of R independent repli-cas of the system is initially prepared at some low (orzero) packing fraction, ϕ . In the simulations reportedbelow R = 10 and ϕ = 0. In parallel, each replicais taken through an annealing schedule, which is a se-quence of increasing packing fractions { ϕ , ϕ , . . . , ϕ K } .Each annealing step ( ϕ i → ϕ i +1 ) consists of two parts.In the first part, an equilibrating procedure is appliedat the initial density, ϕ i . In our study, the equilibratingprocedure is event chain Monte Carlo (ECMC) [9] (seeSec. III B) but other methods such as molecular dynam-ics or Metropolis Monte Carlo would also be suitable.The second part of the annealing step is to increase thepacking fraction from ϕ i to ϕ i +1 holding the relative po-sitions of the particles fixed. Typically some replicas inthe population will now suffer overlaps between spheresand therefore have disallowed configurations. These dis-allowed replicas are removed from the population. Sup-pose that (cid:15)R replicas are culled from the population. Thepopulation size is restored to R by randomly choosing (cid:15)R replicas from among the surviving (1 − (cid:15) ) R replicas andmaking one copy of each of them. For an observable O ,the PA estimator, ˜ O for the ensemble average of an ob-servable, (cid:104)O(cid:105) at each packing fraction is obtained from anaverage over the population at that packing fraction. Thepopulation can be thought of as an approximate, finiterealization of the statistical ensemble for hard spheres.Note that it is possible that there are no allowed config-urations at the new packing fraction. In this case, thealgorithm must be terminated without producing resultsfor higher packing fractions.In the limit of large R , if the original population is anunbiased sample from the equilibrium hard sphere en-semble at packing fraction ϕ i then the new populationwill also be an unbiased sample from the equilibriumhard sphere ensemble at packing fraction ϕ i +1 . How-ever, the new population is now correlated due to thecopying of replicas and, for finite R , the new populationmay be biased due to the omission of configurations thatare important at the higher density but too rare to berepresented in the population at the lower density. The equilibrating procedure at the new density partially cor-rects these problems. An analysis of the error in PA dueto finite population size is given in Sec. II E.Population annealing may be implemented with a fixedannealing schedule however we have found it more conve-nient to use an adaptive annealing schedule in which thefraction culled in each step is fixed. In our simulationsthe culling fraction is (cid:15) = 0 .
05. In this implementation,the annealing schedule is a random list is a random list ofpacking fractions. However, in practice spacing betweensuccessive values of ϕ is very small so it is straightfor-ward to interpolate to obtain observables at any packingfraction. For the adaptive annealing schedule, in princi-ple, the algorithm will never terminate but it may jam inthe sense of taking smaller and smaller steps convergingto a maximal or jammed density.Population annealing for hard spheres in the NVT en-semble is particularly simple because each allowed hardsphere configuration has the same probability in the equi-librium ensemble so the resampling step requires no re-weighting or Boltzmann factors. It would also be possi-ble to use the NPT ensemble (with fixed temperature),in which case a Boltzmann factor exp[ − β ( p i +1 − p i ) V r ]would be required when deciding how many copies tomake of replica r in the step ( p i → p i +1 ) where the vol-ume of replica r at pressure p i is V r . C. Configurational entropy estimator
One desirable feature of PA is that it gives direct ac-cess to thermodynamic potentials. For the hard sphereversion of the algorithm described here, we have directaccess to the configurational entropy S c ( ϕ ) as a functionof packing fraction. From the basic definition of entropy,we have S c = log Ω where Ω is the statistical weight ordimensionless volume in configuration space of accessibleconfigurations and the units are chosen so that Boltz-mann’s constant is unity. The factor 1 − (cid:15) is an estimatorof the ratio of the statistical weight before and after theculling, Ω i +1 / Ω i , during the annealing step ( ϕ i → ϕ i +1 ).Thus we have the following expression for the change inthe estimator ˜ S c in one annealing step,˜ S c ( ϕ i +1 ) − ˜ S c ( ϕ i ) ≈ log(1 − (cid:15) ) . (1)Summing this estimate over the annealing schedule up tostep k , we obtain the entropy at packing fraction ϕ k interms of the entropy at the initial packing fraction andthe number of annealing steps taken,˜ S c ( ϕ k ) = ˜ S c ( ϕ ) + k log(1 − (cid:15) ) . (2)The equation of state can be obtained by differentiatingwith respect to ϕ . These estimates become exact in the R → ∞ limit. Errors are discussed in Sec. II E. D. Weighted Averaging
A very useful feature of PA is the ability to improveboth statistical and systematic errors by combining manyindependent runs using weighted averaging. If the popu-lation is perfectly equilibrated and there are no system-atic errors then the most efficient way to reduce statisticalerrors in an observable is to perform ordinary, unweightedaveraging over the values of the observable in each run.However, for finite R , each run is not completely equili-brated and there are systematic errors. Unweighted aver-aging suppresses statistical errors but does not suppresssystematic errors. However, using weighted averaging wecan also suppress systematic errors. Suppose we havemeasured both the configurational entropy estimator, ˜ S c and the estimator of an observable ˜ O at a given value of ϕ in each of M independent runs, all with the same pop-ulation size. The weighted average for the observable, O is O = (cid:80) Mm =1 ˜ O ( m ) exp[ ˜ S c ( m ) ] (cid:80) Mm =1 exp[ ˜ S c ( m ) ] , (3)where m indexes the independent runs. The observable O must be a quantity that can be measured in a singlereplica at a single packing fraction. Examples includepressure and various correlation functions but not over-laps between replicas or the entropy itself.We make the following claim: Weighted averaging, asdefined in Eq. (3), yields an exact, unbiased result forfixed population size R , and fixed annealing schedule { ϕ , ϕ , . . . , ϕ K } , in the limit of infinitely many runs, M → ∞ . The same conclusion holds for the adaptiveannealing schedule used here and is discussed below.The validity of weighted averaging is trivial for R = 1,(simulated annealing). For a given ϕ i , a run either sur-vives to that packing fraction or is terminated at a lowerpacking fraction. Surviving runs all have the same weightsince there is no resampling, while terminated runs havezero weight. Each surviving run is an unbiased sampletaken from the hard sphere ensemble. Since the weight-ing factor is either 0 or 1, the weighted average is a simpleaverage over the surviving runs. The R = 1 version ofthe algorithm lacks resampling and is therefore highlyinefficient at high densities.The validity of weighted averaging for arbitrary R canbe established using an inductive argument. The popu-lations of the runs at the initial packing fraction ϕ areassumed to be equilibrated so that unweighted averagingis appropriate at ϕ . Correspondingly, the weight fac-tors are the same and the claim is trivially valid for ϕ .Now suppose that the weighted averaging claim holds atpacking fraction ϕ i with weights w ( m ) i . In the first part ofthe resampling step from ϕ i to ϕ i +1 , the random fraction (cid:15) ( m ) i of replicas with hard sphere overlaps are culled fromthe population. Consider the situation before copying isdone to restore the population size to R . Each survivingreplica in run m should still have the weight w ( m ) i if one averages over all replicas in all runs. However, what weactually do is to first average the observable within eachrun and then average over runs. Therefore weight of eachrun must be adjusted to reflect its new population size, w ( m ) i +1 ∝ (1 − (cid:15) ( m ) i ) w ( m ) i , (4)where the constant of proportionality is set by the re-quirement that the weights are normalized. Compar-ing to Eq. (1), we see that if w ( m ) i ∝ exp[ ˜ S c ( m ) ( ϕ i )]then w ( m ) i +1 ∝ exp[ ˜ S c ( m ) ( ϕ i +1 )] and the inductive hypoth-esis is verified. The remainder of the annealing step( ϕ i → ϕ i +1 ) consists of randomly copying replicas torestore the population size to R and then applying theequilibrating procedure to each replica. The copying stepyields additional fluctuations but does not change theexpectation of the observable for each run. The equi-librating procedure may change the expectation of theobservable for a given run if the individual runs were ini-tially out of equilibrium. Nonetheless, the correctness ofweighted averaging is preserved because the equilibratingprocedure and weighted averaging both converge to thesame hard sphere distribution.Weighted averaging for the entropy itself differs fromother observables because the entropy is obtained froma sum over all packing fractions. Nonetheless, the finalresult for the weighted average of the entropy S c takes asimple form, S c = log 1 M M (cid:88) m =1 exp[ ˜ S c ( m ) ] . (5)This result is an example of the Jarzynski equality [22]. Aderivation for PA in the the context of free energy ratherthan the entropy can be found in Ref. [13].Finally, in the adaptive step algorithm used here, theannealing schedule is a random list so each run visits adistinct set of packing fractions whereas it is necessaryto carry out weighted averages at fixed values of ϕ . Toaccomplish this we interpolate both the observables andthe entropy between packing fractions and then use theweighted averaging formulas, Eqns. (3) or (5). E. Equilibration and Systematic Errors
How do we know whether PA is yielding equilibriumvalues of observables? As discussed above, the resam-pling step for finite population size introduces systematicerrors because the distribution at the new packing frac-tion is not fully sampled. The exactness of weighted av-eraging gives a way to quantify the deviations from equi-librium. The analysis of systematic errors is discussed indetail in Ref. [13] and is briefly reviewed here.The expected value of an observable from a single runat population size R of PA is the unweighted averageover runs while the the exact value of the observable isthe weighted average (both averages taken in the limitof infinitely many runs). Thus the systematic error inan observable, ∆ O is given by ∆ O = (cid:104) ˜ O(cid:105) − O where (cid:104) ˜ O(cid:105) represents an unweighted average over independentruns of the algorithm. The difference between weightedand unweighted averages depends on the variance of theweighting factor in Eq. (3). For the case that the jointdistribution of ˜ O and ˜ S c is a bivariate Gaussian, it wasshown in Ref. [14] that systematic errors in measuring O are given by the covariance of ˜ O and ˜ S c ,∆ O = cov( ˜ O , ˜ S c ) = var( ˜ S c ) (cid:34) cov( ˜ O , ˜ S c )var( ˜ S c ) (cid:35) . (6)The second, trivial identity in Eq. 6 is useful becausethe ratio in the square brackets goes to a constant thatdepends on O as R → ∞ while var( ˜ S c ) diminishes as 1 /R and is independent of the observable. The expression onthe far right in Eq. (6) motivates defining an equilibrationpopulation size ρ f , ρ f = lim R →∞ R var( ˜ S c ) . (7)Systematic errors in all observables are proportional to ρ f /R and PA simulations are well-equilibrated when R (cid:29) ρ f . In population annealing ρ f plays the samerole as the exponential autocorrelation time in Markovchain Monte Carlo methods.For sufficiently large R , ˜ S c and other observables arisefrom many independent, additive contributions so non-rigorous “central limit theorem” arguments suggest thatthe joint distribution of entropy and any observableshould be a bi-variate Gaussian. The approach to Gaus-sianity with R was investigated in Ref. [14]. If R istoo small the joint distribution will not be Gaussian andthere will be corrections to Eq. (6) containing higher cu-mulants of the joint distribution. We shall see later thatfor the fixed population size used here, the joint distribu-tion of pressure and entropy is a bivariate Gaussian forsmall packing fractions but significant exponential tailsappear at high packing fractions. These corrections andtheir meaning are explored in Sec. IV AIn our simulations we use weighted averaging exten-sively. What are the systematic errors in weighted av-erages? The same arguments that shows that var( ˜ S c )controls systematic errors in a single run demonstratethat the relevant measure of equilibration for a weightedaverage is var( S c ), where S c is the weighted average esti-mator of the configurational entropy defined in Eq. (5).From var( S c ) we define the weighted average equilibra-tion population size ρ ∗ f ( R ) that is a function of the fixedpopulation size R used in the weighted average, ρ ∗ f ( R ) = lim M →∞ M R var( S c ) . (8)Since weighted averaging using M runs of population size R is less efficient than doing a single large run with pop-ulation size M R , we have that ρ ∗ f ( R ) ≥ ρ f and we alsoexpect that lim R →∞ ρ ∗ f ( R ) → ρ f . Of course, it is not practical to compute var( S c ) fromrepeated experiments each with M runs. Instead we usea bootstraps procedure to estimate var( S c ) by resamplingwith replacement from our M runs. Each sample of size M is employed to compute S c from Eq. (5) and the vari-ance is compute by resampling many times.It is known that in the thermodynamic limit and atsufficiently high packing fraction the equilibrium state ofthe binary mixture studied here is phase separated withseparate crystals of the small and large spheres [23]. Fora finite number of spheres the dominant configurationsin the equilibrium distribution are not known since theentropic cost of the interface between the two crystallineregions may prevent phase separation. It may be that theequilibrium states are dominated either by phase sepa-rated configurations or other non-random configurations.We have not seen evidence of either phase separation orother forms of ordering in our PA simulations. F. Statistical Errors
If an observable O is averaged over R independentmeasurements, the statistical error δ O in the mean isgiven by (cid:112) var( O ) /R . The resampling step in PA in-troduces correlations so that statistical errors are largerthan (cid:112) var( O ) /R . As discussed in detail in Ref. [14], wecan bound statistical errors by ignoring the de-correlatingeffects of the equilibrating procedure. If there is no equili-brating procedure, each descendent of an initial replica isthe same. We refer to the set of descendants of an initialreplica as a family and define n i to be the fraction of thepopulation in family i . In the absence of the equilibrat-ing procedure O takes a single value for every member ofthe family, call this value O i . Given this assumption,˜ O = (cid:88) i O i n i . (9)To proceed we make two additional assumptions: (1) theobservable O i and family fraction n i are uncorrelated and(2) the distribution of family fraction n i has small fluctu-ations from run to run for large R . Given these assump-tions, the variance of the observable is bounded by,var( ˜ O ) ≤ var( O ) (cid:88) i n i . (10)Assuming that the second moment of the family frac-tion scales as 1 /R we define the mean square family size , ρ t , ρ t = lim R →∞ R (cid:88) i n i . (11)The bound on the statistical error in δ ˜ O becomes δ ˜ O ≤ (cid:112) var( O ) ρ t /R. (12)Note that if there is no decimation so that n i = 1 /R ,then ρ t = 1 and Eq. (12) (as an equality) reduces tothe expression for statistical errors for uncorrelated mea-surements. Since ρ t provides only a bound on statisticalerrors, in practice we use bootstrapping to estimate er-rors.In Ref. [24] it is shown that ρ f and ρ t are close toone another. We investigate the relation between thesetwo measures in Sec. IV A. Since measuring ρ f requiresmultiple runs while ρ t can be measured from a single run, ρ t can serve as a practical measure of equilibration.It is important to note that in the regime that R (cid:29) ρ f ,systematic errors diminish as 1 /R while statistical errorsdiminish as 1 / √ R so that, generally, systematic errorsare much smaller than statistical errors. G. Comparison to Parallel Tempering
Population annealing is a sequential Monte Carlomethod [25] in contrast to most simulation methodsin statistical physics, which are Markov chain MonteCarlo methods. In sequential Monte Carlo, equilibra-tion is approached as the population size R increaseswhile for Markov Chain Monte Carlo, equilibration isapproached as the number of Monte Carlo sweeps is in-creased. Among Markov chain Monte Carlo methods par-allel tempering shares the greatest similarity with pop-ulation annealing. It is comparably efficient to paralleltempering but, as we have seen, it has some interestingadvantages–it gives direct access to thermodynamic po-tentials such as entropy or free energy, it is amenableto massive parallelization and multiple runs can be com-bined to improve both statistical and systematic errors.Parallel tempering explores the disconnected minima ina rough free energy landscapes by annealing to low tem-peratures multiple times, using correctly designed swapmoves to insure that each region of the free energy land-scape is visited with the statistical weight determined bythe Gibbs distribution. Population annealing exploresrough free energy landscapes in a single annealing runusing a large number of replicas to populate the discon-nected minima in the landscape. Resampling insures thateach minimum contains a fraction of the population givenby the Gibbs distribution.For parallel tempering and other Markov chain MonteCarlo methods, systematic and statistical errors canbe determined from autocorrelation functions of observ-ables. Systematic errors are related to the exponentialautocorrelation time while statistical errors are relatedto the integrated autocorrelation time. Thus the twocharacteristic population sizes, ρ f and ρ t , play analogousroles for PA to the exponential and integrated autocor-relation times, respectively. III. MODEL AND METHODSA. Hard Sphere Model
Hard spheres have long offered a simple and use-ful model for studying the properties of fluids. Al-though these systems have been studied for decades,there are still mysteries surrounding some of their behav-iors. Monodisperse systems exhibit a distinct entropy-driven first-order phase transition from the disorderedliquid to a crystalline solid, but a metastable fluid branchpersists beyond the transition and its high-density be-havior is still not fully understood. It is as yet unclearto what packing fraction this metastable branch per-sists [7, 11, 18, 19], and there has also been evidencefor [18, 26] and against [8, 11] the existence of a thermo-dynamic glass transition.One difficulty in studying the high-density behaviorof the metastable fluid is that, since the true thermo-dynamic state is a solid, full equilibration will lead tocrystallization. Size polydispersity is often introduced toavoid crystallization. We study 50/50 mixture of sphereswith diameter ratio 1.4:1, which does not easily crystal-lize [18, 23, 27].At low and moderate densities the equation of stateof a binary mixture of hard spheres is accuratelydescribed by the Boubl´ık-Mansoori-Carnahan-Starling-Leland (BMCSL) equation[20, 21]: Z ( ϕ ) = (1 + ϕ + ϕ ) − ϕ ( y + y ϕ ) − y ϕ (1 − ϕ ) (13)The constants y , y , and y depend on the choice ofpolydispersity and for the 50/50 binary mixture withpolydispersity ratio 1.4:1 they are given by, y = 0 . y = 0 . y = 0 . B. Event Chain Monte Carlo
Population annealing requires an equilibrating proce-dure performed at each annealing step. We use eventchain Monte Carlo (ECMC), which has been shown tobe highly efficient for simulating 2D and 3D hard spheresystems [9, 10, 28–31].In an ECMC step, a particle is chosen at random andmoved in a given direction. The particle moves until itcollides with some other particle (the “event”). The orig-inal particle remains at the point of collision, while theparticle that was struck then moves in the same direc-tion. This process continues until the total displacementlength of the chain of particles reaches a predeterminedchain length (cid:96) c . The chain length is a parameter thatcan be adjusted, e.g. to minimize correlation times.There are several variants of ECMC, but the simplestand fastest [9] version has only two directions of mo-tion, + x or + y . This method, called “x-y straight eventchain,” violates detailed balance, but preserves globalbalance and thus converges to the correct equilibriumstate [9]. C. Population Annealing Pseudo-Code
The following is pseudocode for our implementation ofpopulation annealing with event chain Monte Carlo fora polydisperse hard sphere system. The diameters of thesmall and large particles are fixed respectively at 1 and1.4.1. Initialize each of R replicas at zero packing fractionby choosing N points at random within a unit cubesimulation cell with periodic boundary conditions.Each particle is also labeled with a size accordingthe prescription for polydispersity.2. For each replica, determine d min , the minimumrelative center-to-center distance between any twoparticles as shown in Fig. 1: d min = min i,j (cid:20) | r i − r j | σ i + σ j (cid:21) , (14)where r i is the position of particle i and σ i is itsdiameter. The dimensions of this system and par-ticle positions may be rescaled by 1 /d min withoutcausing overlaps. Note that for the first annealingstep this rescaling increases the size of the system,and for all later steps the size is decreased.3. Sort the list of d min values and define the ( (cid:15)R ) th smallest value of d min as d ∗ .4. Rescale the simulation box size L , packing fraction ϕ and position r of each particle according to L ← L/d ∗ , ϕ ← ϕ/ ( d ∗ ) and r ← r /d ∗ . (Note that ϕ isinitially zero and needs to be calculated after thefirst annealing step).5. After rescaling there are (cid:15)R replicas with invalidconfigurations [32]. These replicas are eliminatedfrom the population. The same number of repli-cas are randomly chosen (with replacement) fromamong the valid replicas and copies of these repli-cas are added to the population. After this step,the population consists of R valid configurations.6. Perform one ECMC sweep on each replica. A sweepshould displace on average N particles. We choosethe chain length (cid:96) c to displace an average of √ N particles (See Appendix A) and then perform √ N ECMC moves.7. Repeat steps 2-6 on each replica until a predefinedpacking fraction is reached.
FIG. 1: The minimum relative center-to-center distance overall particle pairs in a replica is d min . D. Observables
1. Entropy
The relative configurational entropy ˜ S c at each stepin the annealing schedule is obtained from Eq. (2). Thesimulations are initialized at ϕ = 0 and the entropyat zero packing fraction is taken to be zero so that allentropy values are negative. We would like to know theweighted average value of the entropy ¯ S c as a function of ϕ . However, our annealing schedule traverses a fixed setof entropies and the packing fraction at each entropy isvariable. Values of the entropy at each packing fractionare obtained by interpolation and Eq. (5) is applied toobtain the weighted average entropy.
2. Equation of State
The dimensionless equation of state, Z is expressed interms of the pressure P , Z = πσ Pϕ , (15)where σ is the average cubed sphere diameter, which fora binary mixture with diameters σ A and σ B is given by σ = ( σ A + σ B ). We obtain Z at each annealing stepand then use interpolation to determine Z at a fixed set ofpacking fractions. We report the weighted average equa-tion of state ¯ Z as a function of packing fraction. Errorbars are obtained by bootstrapping the weighted aver-ages. Two independent ways of measuring the pressureare described in the next two subsections.
3. Thermodynamic Pressure Estimator
Using Eq. (2), pressure and equation of state can befound through the thermodynamic definition, PT = ∂S∂V . (16)In our simulations the change in entropy for each an-nealing step is a constant log(1 − (cid:15) ) and the incrementin packing fraction is variable. Thus a thermodynamicestimator for the pressure at annealing step k is P k T = − ϕ k V k ∆ S ∆ ϕ k = − ϕ k V k log(1 − (cid:15) )∆ ϕ k . (17)Setting T = 1 , Z = − πσ V k log(1 − (cid:15) )∆ ϕ k (18)where V k is the volume at step k and ∆ ϕ k = ( ϕ k +1 − ϕ k − ) /
4. Dynamic Pressure Estimator
Event chain Monte Carlo offers a direct method of cal-culating the equation of state, as described in Ref. [29].Consider an event chain in the + x direction. At eachcollision between particles j and k , the distance betweenthe centers of the particles projected on the direction ofmotion is x k − x j . Define the “lifted” distance of an eventchain, x final − x initial as x final − x initial = (cid:96) c + (cid:88) ( k,j ) ( x k − x j ) (19)where the sum is over all of the collision events in thechain. The equation of state is then given by an averageof the lifted distance over all event chains [29], Z = (cid:28) x final − x initial (cid:96) c (cid:29) chains . (20)The results we present below for the equation of state arefrom this dynamic, ECMC estimator.We note that Eq. (20) is not a correct measure of Z when (cid:96) c > L . When using the dynamic event chainlength described in the Appendix, Eq. (20) cannot beused for ϕ (cid:46) . E. ρ f , ρ ∗ f and ρ t In order to assess equilibration we computed the pop-ulation size measures ρ f , ρ ∗ f , and ρ t as described in Secs.II E and II F. The mean square family size ρ t is straight-forward to measure in each run with Eq. (11) by keepingtrack of the family to which each replica belongs. Theseestimates can also be combined through a weighted av-erage across independent runs.To estimate the equilibrium population size ρ f throughEq. (7), however, we need the variance of entropy overmany independent runs. However, because we use a vari-able annealing schedule, we do not have direct access tothe variance of the entropy as a function of packing frac-tion. Instead we have variable packing fractions at eachvalue of the entropy and var( ϕ ) as a function of S c . By considering the S c vs. ϕ curves from each run as inde-pendent random functions we can make the first-orderestimate var( ˜ S c ) = (cid:18) ∂S c ∂ϕ (cid:19) var( ϕ ) , (21)with the same symmetric derivative estimator for ∂S c /∂ϕ as in Eq. (17).The situation is simpler for the weighted average equi-librium population size ρ ∗ f , because the weighted averageentropy ¯ S c is calculated at fixed packing fractions, sovar( ¯ S c ) can be calculated directly by bootstrapping andEq. (8) gives ρ ∗ f . IV. RESULTS
In this section we present the results of our simulations.Our results address two objectives. The first objective isto study the properties of the population annealing al-gorithm as applied, for the first time, to fluid systemsand determine whether the algorithm can produce equi-librated results in the glassy regime of bidisperse hardspheres. The second goal is to obtain precise equilibratedresults for the equation of state at high densities andto study nonequilibrium behavior at even higher densi-ties near jamming for the 50/50 binary mixture of hardspheres with diameter ratio 1.4:1. We study systems with N = 60 and N = 100 particles. For each of these caseswe performed M ≈ R = 10 . Table I shows the simulationparameters. Note that the population size (10 ) and thenumber ECMC sweeps per annealing step (one) are thesame for both N = 60 and 100 particles but the numberof annealing steps required to keep the culling fractionfixed increases linearly in the number of particles. Thislinear scaling is the result of the fact that the probabilityof an overlap for a given fractional compression increas-ing linearly in the increase in packing fraction. Togetherwith the linear scaling of carrying out a single sweep ofECMC, we see that the naive complexity of the algorithmis O ( N ). Of course, this scaling does not take into ac-count how computational resources must scale with N to achieve equilibration. The question of equilibration isdiscussed in the following subsection. A. Equilibration of Population Annealing
Figure 2 shows the characteristic population sizes ρ f and ρ t vs. packing fraction. These quantities are dis-cussed in Secs. II E and II F, respectively, and they char-acterize the errors in population annealing. Specifically,when ρ f is small compared to the population size, sys-tematic errors are small and each run is well-equilibrated.Similarly, when ρ t is small compared with the populationsize, statistical error are small. Figure 2 demonstrates TABLE I: Simulation parameters for the results reported inSec. IV. ϕ max is the maximum packing fraction. Average wallclock time for a run corresponds to an OpenMP implementa-tion with 32 threads.System Size N
60 100Population Size R Max packing fraction ϕ max (cid:15) M ρ f (solid lines) and ρ t (dotted lines) vs. packing fraction ϕ forsizes N = 60 (red) and N = 100 (blue). See Secs. II E andII F for definitions. that both measures are relatively flat and much smallerthan the population size ( R = 10 ) until ϕ = 0 .
55, butthen both ρ f and ρ t grow quickly, increasing by abouta factor of 100 between ϕ = 0 .
55 and 0.58. If we im-pose a relatively conservative equilibration criterion that
R > ρ f then individual runs are found to be well-equilibrated for ϕ < .
575 for both system sizes.Note that ρ t is approximately a factor of 2 larger than ρ f at low density but in the region 0 . (cid:46) ϕ (cid:46) .
59, thetwo quantities track each other closely (see Ref. [24]). Fi-nally, for ϕ ≈ .
59 we find that ρ f ≈ ρ t ≈ R . Beyondthis packing fraction, individual runs are clearly not inequilibrium. Since ρ t is estimated from runs at popula-tion size R , it bounded by R while ρ f is not bounded sothe two quantities part ways above ϕ (cid:38) . ρ f for the two system sizes itseems clear that in the glassy regime ϕ ≥ .
58 it is moredifficult to equilibrate larger systems. With only two
TABLE II: Measured quantities at selected packing fractions, ϕ for N = 60 spheres. ¯ S c is the weighted average entropy de-fined in Eq. (5) and ¯ Z is the weighted average equation ofstate with one standard deviation errors. µ , σ , and λ arethe parameters of the best fit exponentially modified Gaus-sian distribution (see text above Eq. (23)). A dash in the λ column indicates that a Gaussian distribution is the preferredfit. ρ f and ρ ∗ f are the equilibration population sizes for sin-gle runs and weighted averages defined in Eqs. (7) and (8),respectively. ϕ ¯ S c ¯ Z µ σ λ ρ f ρ ∗ f × × × × × × × × × × TABLE III: Measured quantities at selected packing frac-tions, ϕ for N = 100 spheres. See Table II for details. ϕ ¯ S c ¯ Z µ σ λ ρ f ρ ∗ f × × × × × × × × × × system sizes, we cannot determine the scaling behaviorof ρ f with N .Results from weighted averaging are equilibrated tohigher densities than results from a single run. Figure 3shows ρ ∗ f ( R ) the weighted average equilibration popula-tion size, together with ρ f , as a function packing fraction.If we apply the conservative equilibration criterion that M R ≥ ρ ∗ f ( R ) then results from weighted averagingare well-equilibrated for ϕ < .
587 for N = 100 hardspheres and for ϕ < .
591 for N = 60 hard spheres.Note that we believe that our results for the fluid equa-tion of state remain reasonably accurate to somewhathigher packing fractions than these strict cut-offs.A deeper understanding of weighted averaging can begained by looking at the joint distribution of the entropyand equation of state estimators, ˜ S c and ˜ Z , as a functionof packing fraction. Figure 4 displays these joint distri-butions as scatter plots for N = 60 particles and severalpacking fractions. Each point in these figures representsa single run of PA. The horizontal coordinate of the pointis the relative configurational entropy estimator ˜ S c andthe vertical coordinate is the equation of state estima-tor ˜ Z for the run. The x and y coordinates of the redsquares are the weighted average values ¯ S c and ¯ Z , re-spectively. Two important features of these plots are im-mediately evident: (1) there is an approximately inversecorrelation between ˜ S c and ˜ Z so that the weighted aver-age value of the pressure is less than the ordinary averageand (2) there are many outliers with large entropies and FIG. 3: (color online) Equilibration population sizes, ρ f (solidlines) and ρ ∗ f (dotted lines) vs. packing fraction ϕ for sizes N = 60 (red) and N = 100 (blue). See Sec. II E for definitions.Dotted horizontal lines at 10 and 10 show where the cut-offsfor equilibration occur for single runs and weighted averages,respectively, based on the requirements 0 . R > ρ f for singleruns and 0 . MR > ρ ∗ f for weighted averages. small pressures and the distributions are clearly not bi-variate Gaussians. For ϕ = 0 .
56 and 0.58, these outliershave little effect on the weighted averages but for the ϕ = 0 . S c distributions (bluecurves) for ϕ = 0 .
58, and 0.60. The distributions havebeen smoothed with Gaussian kernels using Mathemat-ica’s
KernelMixtureDistribution . It is clear that these dis-tributions have exponential tails for large entropy ratherthan Gaussian tails. A good fit to these distributions,shown as dotted curves in Fig. 5 , is obtained using anexponentially modified Gaussian, which is defined as fol-low: Let X be Gaussian distributed and Y be exponen-tially distributed, then Z = X + Y is described by anexponentially modified Gaussian distribution. An expo-nentially modified Gaussian requires three parameters: µ and σ are the mean and standard deviation of the Gaus-sian distribution, respectively, and λ is the characteristicdecay rate of the exponential distribution. The probabil-ity density function is given by p ( x ; µ, σ, λ ) = λ e λ ( λσ − x +2 µ )erfc (cid:18) λσ − x + µ √ σ (cid:19) . (22)The best fit values of the parameters for several valuesof ϕ are shown in Tables II and III, for N = 60 and 100particles, respectively. FIG. 4: (color online) Joint distribution of the entropy andequation of state estimators for N = 60 and packing fractions ϕ = 0.56, 0.58, 0.60, and 0.62 for Figs. 4a-d, respectively.Each point represents a single run of population annealing.The x -coordinate of the point is the entropy estimator, ˜ S c and the y -coordinate is the equation of state estimator, ˜ Z .The weighted average values are shown as red squares.FIG. 5: Distributions of the configurational entropy esti-mator ˜ S c (solid curves) at packing fractions ϕ = 0 .
58 (Fig.5a) and ϕ = 0 .
60 (Fig. 5b) together with a fit to an expo-nentially modified Gaussian distribution (dotted curves) for N = 60 particles. See Table II for fitting parameters. S c , defined in Eq. (5), approaches, as M → ∞ , the inte-gral of the distribution and the weighting factor,¯ S c = log (cid:90) dx p ( x ; µ, σ, λ ) e x . (23)Carrying out the integral we find,¯ S c → µ + σ (cid:18) λλ − (cid:19) . (24)When the exponential decay parameter λ is significantlylarger than one, the effect of the exponential tail is smalland the weighted average is a good approximation for theentropy and other observables so long as the number oftrials is sufficiently large to explore the tail, as is the casehere with M ≈ λ (cid:46) λ (cid:46)
1, the population annealing results do notyield useful information about the equilibrium propertiesof the system. From Tables II and III we see that for N =60, λ = 1 is reached at about ϕ = 0 .
60 while for N = 100, λ = 1 is reached between 0.59 and 0.60. These estimatesfor where equilibration breaks down are consistent withthe more conservative estimates based on ρ ∗ f . B. Equation of State
Figure 6 shows the weighted average equation of state¯ Z as a function of packing fraction ϕ for systems with N = 60 and N = 100 particles. The weighted average isperformed using M ≈ R = 10 . We estimated thestatistical errors by randomly resampling with replace-ment (bootstrapping) the collection of trials and recalcu-lating the weighted average for the resampled trials. Wenote that this method makes use of the joint distributionof entropy and pressure to estimate systematic errors inthe equation of state. The shaded regions around theweighted average curves represent 95% confidence inter-vals. If the original weighted average is dominated bya small number of trials with significantly greater en-tropies, the bootstrapped collection may leave out thesetrials, resulting in a highly skewed confidence interval.This effect is clearly seen in Fig. 6 at ϕ ≈ .
61 where thesimulations have clearly fallen out of equilibrium.The solid line in Fig. 6 is the BMCSL equation of state(see Eq. (13)), which is reasonably accurate for low densi-ties but significantly underestimates Z for high densities.Of course, the BMCSL equation predicts finite pressureup to ϕ = 1 so that deviations of this kind are inevitablebut it is not clear where they become significant. Figure FIG. 6: (color online) The weighted average equation of state,¯ Z with bootstrapped 95% confidence intervals for N = 60(red) and N = 100 (blue) particles. The solid black line isthe BMCSL equation of state.FIG. 7: (color online) Deviations ∆ Z between the simulationresults and the BMCSL equation of state as a function ofpacking fraction ϕ with bootstrapped 95% confidence inter-vals for N = 60 (red) and N = 100 (blue) particles. Z = ¯ Z − Z CS between the equa-tion of state measured in the simulations and the BM-CSL prediction. For N = 60, significant deviations beginat ϕ = 0 .
56 while for N = 100 small deviations persistto much lower packing fractions. Large differences fromBMCSL begin at ϕ = 0 . g ( r ). The pair correlation function is shown in Fig. 8 for N = 100 particles at packing fraction ϕ = 0 .
58. Thepair correlation function for N = 60 particles is indis-tinguishable from g ( r ) for N = 100 particles. We seesharp peaks only at the three contact distances for thisbinary mixture. If the system typically formed phase-separated crystals or other ordered states we would ex-pect additional sharp structure in g ( r ). In addition, Fig.8 is obtained from a weighted average over runs but weobserved no noticeable difference between the weighted1 FIG. 8: Pair correlation function g ( r ) for N = 100 particlesat packing fraction ϕ = 0 .
58 as a function of separation r .The three peaks correspond to the separations at contact forthis binary mixture.FIG. 9: (color online) Relative difference δZ/Z between thedynamic and thermodynamic estimates of the equation ofstate, Eqs (20) and (17), respectively, for N = 60 (red) and N = 100 (blue) particles. and unweighted averages. This result suggests that theentropy of these systems is not strongly correlated withthe g ( r ). Although we cannot rule out the possibility thata small fraction of configurations in the ensemble are or-dered, we conclude that the predominant configurationssampled by the algorithm are amorphous.The above results for the equation of state have beenobtained using the dynamic measure of the pressure.Figure 9, shows the relative discrepancy, δZ/Z betweenthe thermodynamic and dynamic measures of pressure,given in Eqs. (17) and (20), respectively. We find that δZ/Z ≈ − until ϕ ≈ .
59. The good agreement be-tween these two independent measures is an importantvalidation of the algorithm and gives us confidence thatestimates of the equation of state based on Eq. (17) arequite accurate.
TABLE IV: Parameters and bootstrapped 95% confidenceinterval for fitting the equation of state to the free volumeform of Eq. (25) for the range of 0 . < ϕ < . N d (cid:48) ϕ c
60 2 . +0 . − . . +0 . − .
100 2 . +0 . − . . +0 . − . FIG. 10: High density non-equilibrium equation of state, ¯ Z (solid curves) along with fits (dotted curves) to the free vol-ume form, Eq. (25) for N = 60 (red) and N = 100 (blue)particles. C. High Density Behavior
Though the simulations have fallen out of equilibriumfor ϕ (cid:38) .
6, population annealing continues to work as anonequilibrium protocol to achieve high density packingsand nearly jammed states. We fit our equation of statedata in the range 0 . < ϕ < .
625 to the “free volume”form, Z = d (cid:48) ϕ c ϕ c − ϕ . (25)Table IV gives the parameters of the fit and Fig. 10shows the simulation data along with fits. It should benoted that in this range of packing fractions the weightedaverage is dominated by the single highest entropy run sothat it is likely that the estimate for ϕ c would increaseas the number of runs or population size increases. Itwould be interesting to study how ϕ c changes as eitherpopulation size or number of runs increases. Our resultsagree reasonably well with Odriozola and Berthier [8]who find d (cid:48) = 2 .
82 and ϕ c = 0 . . < ϕ (cid:46) .
65 .To explore the behavior of population annealing as aprotocol for achieving nearly jammed states, we allowedten runs to go to much higher packing fraction. Theseruns terminated after exceeding a maximum run timeand reached packing fractions in the range 0 . < ϕ c < . ϕ c in the range 0 . < ϕ c < . ϕ c are smaller than the estimates shown in Ta-ble IV because they represent only 10 rather than 1000independent experiments. V. DISCUSSION
We have developed population annealing Monte Carlofor fluid systems and applied it to a glass-forming binarymixture of hard spheres. We find that population an-nealing is a promising method for computational studiesof fluids in the high density regime. It is a highly par-allelized algorithm that is well-suited to simulations onlarge computer clusters either by using weighted averag-ing and a large number of independent runs, as we havedone here, or by carrying out a smaller number of verylarge population runs with a massively parallelized imple-mentation of the algorithm. The advantage of the formerapproach is that we can learn much about the equilibra-tion of weighted averages from the statistics of multipleruns but this advantage is balanced by the lesser effi-ciency of weighted averaging as opposed to a single largepopulation run.We measured the equation of state in the glassy regionof a 50/50 mixture of hard spheres with diameter ratio1.4:1 and have obtained precise results that are in reason-able agreement with previous simulations using paralleltempering [8]. We find good agreement with the BM-CSL equation of state up to a packing fraction of 0.58but strong deviations above that packing fraction. Al-though this is also the region where equilibration becomesmuch for difficult, we believe the simulations are reason-ably well equilibrated up to a packing fraction of 0.60.We also studied the equation of state at higher densitieswhere population annealing serves as a non-equilibriumprotocol for generating nearly jammed states.We must emphasize that our results at high densitieslikely describe the metastable equilibrium fluid branchof the equation of state and not the true equilibrium,which presumably consists either of a phase separatedcrystalline state or, for small N , perhaps a non-randombest packing of N spheres. The definition of metastableequilibrium depends on the protocol used to sample thedistribution so, in the region where the true equilibriumensemble contains a significant contribution from non-random states, our results may differ from those obtainedfrom other algorithms. It would be interesting to inves-tigate the true equilibrium states of the binary mixturestudied here for finite N .This work is the first application of population anneal-ing to classical fluids and there are undoubtably manyavenues for improvement. We have implemented the al-gorithm in the NVT ensemble but it is important to studypopulation annealing in the NPT ensemble. It would alsobe interesting to explore other annealing schedules. Forexample, a variable number of Monte Carlo sweeps withmore sweeps in the region of the glass transition mayimprove efficiency. Appendix A: Dynamic chain length
Based on the average number of particles displaced ineach event chain, we define an ECMC sweep to be a num- ber of event chain moves that will move approximately N particles on average. We arbitrarily choose to divide thiscomputational work so that an average of √ N particlesare displaced in each event chain. Then √ N event chainsconstitutes an ECMC sweep of the system.The event chain length, (cid:96) c must be chosen dynamically,which can be done using Eq. (20), re-written in the form, Z = 1 + (cid:10) (cid:80) chains ( x k − x j ) (cid:11) chains (cid:96) c . (A1)The “lifting distance” x j − x i between two particles ofdiameter σ with bond orientation ( θ, ϕ ) at contact is x j − x i = σ sin θ cos ϕ. (A2)In a fluid, this bond orientation should be random, sothe average projected distance is (cid:10) x j − x i (cid:11) = 12 π (cid:90) π (cid:90) π/ − π/ d Ω σ sin θ cos ϕ = σ . (A3)This will be somewhat less accurate in an fcc solid, butit turns out to be close enough to give a reasonable esti-mate. If an event chain consists of n c collisions, (cid:10) (cid:88) chains ( x j − x i ) (cid:11) chains ≈ n c σ/ , (A4)and we find, Z ≈ (cid:96) c + n c σ/ (cid:96) c . (A5)We would like ¯ n c ≈ √ N . For density ϕ i then we usethe equation of state Z ( ϕ i − ) at the previous density toestimate (cid:96) c ( ϕ i ) as, (cid:96) c ( ϕ i ) ≈ σ √ N / Z ( ϕ i − ) − . (A6)Empirically, we found a somewhat better approximationto n c = √ N using (cid:96) c ( ϕ i ) = σ (cid:112) N/ Z ( ϕ i − ) − . (A7) Acknowledgments
This work was supported by the National ScienceFoundation (Grant No. DMR-1507506). We thank theMassachusetts Green High Performance Computing Cen-ter (MGHPCC) and the University of MassachusettsAmherst for providing computing resources.3 [1] R. H. Swendsen and J.-S. Wang,
Replica Monte Carlosimulations of spin glasses , Phys. Rev. Lett. , 2607(1986).[2] K. Hukushima and K. Nemoto, Exchange Monte Carlomethod and application to spin glass simulations , J. Phys.Soc. Jpn. , 1604 (1996).[3] H. G. Katzgraber, M. Palassini, and A. P. Young, MonteCarlo simulations of spin glasses at low temperatures ,Phys. Rev. B , 184422 (2001).[4] R. A. Banos, A. Cruz, L. A. Fernandez, J. M. Gil-Narvion, A. Gordillo-Guerrero, M. Guidetti, A. Maio-rano, F. Mantovani, E. Marinari, V. Martin-Mayor, et al., Nature of the spin-glass phase at experimental lengthscales , Journal of Statistical Mechanics: Theory and Ex-periment , P06026 (2010).[5] B. Yucesoy, H. G. Katzgraber, and J. Machta,
Evi-dence of non-mean-field-like low-temperature behavior inthe Edwards-Anderson spin-glass model , Phys. Rev. Lett. , 177204 (2012).[6] T. Okabe, M. Kawata, Y. Okamoto, and M. Mikami,
Replica-exchange Monte Carlo method for the isobaric-isothermal ensemble , Chemical Physics Letters , 435(2001).[7] G. Odriozola,
Replica exchange Monte Carlo applied tohard spheres , J. Chem. Phys. , 144107 (2009).[8] G. Odriozola and L. Berthier,
Equilibrium equation ofstate of a hard sphere binary mixture at very large den-sities using replica exchange Monte Carlo simulations , J.Chem. Phys. , 054504 (2011).[9] E. P. Bernard, W. Krauth, and D. B. Wilson,
Event-chain Monte Carlo algorithms for hard-sphere systems ,Phys. Rev. E , 5 (2009).[10] M. Isobe and W. Krauth, Hard-sphere melting and crys-tallization with event-chain Monte Carlo , J. Chem. Phys. , 084509 (2015).[11] L. Berthier, D. Coslovich, A. Ninarello, and M. Ozawa,
Equilibrium sampling of hard spheres up to the jammingdensity and beyond , Phys. Rev. Lett. , 238002 (2016).[12] K. Hukushima and Y. Iba, in
The Monte Carlo Method InThe Physical Sciences: Celebrating the 50th Anniversaryof the Metropolis Algorithm , edited by J. E. Gubernatis(AIP, 2003), vol. 690, pp. 200–206.[13] J. Machta,
Population annealing with weighted averages:A Monte Carlo method for rough free-energy landscapes ,Phys. Rev. E , 026704 (2010).[14] W. Wang, J. Machta, and H. Katzgraber, Population an-nealing: Theory and application in spin glasses , Phys.Rev. E , 063307 (2015).[15] W. Wang, J. Machta, and H. Katzgraber, Evidenceagainst a mean-field description of short-range spinglasses revealed through thermal boundary conditions ,Phys. Rev. B , 184412 (2014).[16] W. Wang, J. Machta, and H. G. Katzgraber, Chaosin spin glasses revealed through thermal boundary con-ditions , Phys. Rev. B (2015). [17] W. Wang, J. Machta, and H. G. Katzgraber,
ComparingMonte Carlo methods for finding ground states of Isingspin glasses: Population annealing, simulated annealing,and parallel tempering , Phys. Rev. E , 013303 (2015).[18] L. Berthier and T. Witten, Glass transition of dense flu-ids of hard and compressible spheres , Phys. Rev. E ,021502 (2009).[19] G. Brambilla, D. El Masri, M. Pierno, L. Berthier,L. Cipelletti, G. Petekidis, and A. B. Schofield, Probingthe equilibrium dynamics of colloidal hard spheres abovethe mode-coupling glass transition , Phys. Rev. Lett. ,085703 (2009).[20] T. Boublk,
Hard-sphere equation of state , J. Chem. Phys. , 471 (1970).[21] G. A. Mansoori, N. F. Carnahan, K. E. Starling, andT. W. Leland, Equilibrium thermodynamic properties ofthe mixture of hard spheres , J. Chem. Phys. , 1523(1971).[22] C. Jarzynski, Equilibrium free-energy differences fromnonequilibrium measurements: A master-equation ap-proach , Phys. Rev. E , 5018 (1997).[23] A. B. Hopkins, F. H. Stillinger, and S. Torquato, Densestbinary sphere packings , Phys. Rev. E , 021130 (2012).[24] C. Amey and J. Machta, Optimized population annealingfor spin glasses , in preparation.[25] A. Doucet, N. de Freitas, and N. Gordon, eds.,
SequentialMonte Carlo Methods in Practice (Springer-Verlag, NewYork, 2001).[26] M. Hermes and M. Dijkstra,
Thermodynamic signatureof the dynamic glass transition in hard spheres , J. Phys.:Consens. Matter , 104114 (2010).[27] C. O’Hern, S. A. Langer, A. J. Liu, and S. R. Nagel, Ran-dom packings of frictionless particles , Phys. Rev. Lett. , 075507 (2002).[28] M. Engel, J. A. Anderson, S. C. Glotzer, M. Isobe, E. P.Bernard, and W. Krauth, Hard-disk equation of state:First-order liquid-hexatic transition in two dimensionswith three simulation methods , Phys. Rev. E , 042134(2013).[29] M. Michel, S. C. Kapfer, and W. Krauth, General-ized event-chain Monte Carlo: Constructing rejection-free global-balance algorithms from infinitesimal steps , J.Chem. Phys. , 054116 (2014).[30] E. Bernard, Ph.D. thesis, Universite Pierre et MarieCurie (2011).[31] E. P. Bernard and W. Krauth,
Two-step melting in twodimensions : First-order liquid-hexatic transition , Phys.Rev. Lett. , 155704 (2011).[32] Note that this choice of rescaling introduces a tiny de-viation from the hard sphere distribution because onepair of particles out in one replica is alwaysalways