Monte Carlo methods for estimating depletion potentials in highly size-asymmetrical hard sphere mixtures
MMonte Carlo methods for estimating depletion potentials in highlysize-asymmetrical hard sphere mixtures
D. J. Ashton, V. S´anchez-Gil, and N. B. Wilding Department of Physics, University of Bath, Bath BA2 7AY, U.K. Instituto de Qu´ımica F´ısica Rocasolano, CSIC, Serrano 119, E-28006 Madrid,Spain
We investigate Monte Carlo simulation strategies for determining the effective (“depletion”) potential betweena pair of hard spheres immersed in a dense sea of much smaller hard spheres. Two routes to the depletionpotential are considered. The first is based on estimates of the insertion probability of one big sphere in thepresence of the other; we describe and compare three such methods. The second route exploits collective(cluster) updating to sample the depletion potential as a function of the separation of the big particles; wedescribe two such methods. For both routes we find that the sampling efficiency at high densities of smallparticles can be enhanced considerably by exploiting ‘geometrical shortcuts’ that focus the computationaleffort on a subset of small particles. All the methods we describe are readily extendable to particles interactingvia arbitrary potentials.
I. INTRODUCTION
Effective potentials arise in theories of complex multi-component fluids such as colloidal suspensions or poly-mer solutions which comprise mixtures of big and smallparticles. For such a system one seeks to integrate outfrom the full Hamiltonian the degrees of freedom of thesmall particles in order to obtain an effective Hamilto-nian for the big particles. The motivation for doing sois to create an approximate, yet analytically tractabledescription of the true system in terms of a single com-ponent model of big particles. Unfortunately, obtainingthe full effective Hamiltonian is a tall order . A firststep in any theoretical treatment is therefore to deter-mine the two-body effective potential between a singlepair of the big particles in a sea of the smaller species.However, even this task is challenging when there existsa very large disparity in size between the particles asis common in suspensions containing a mixture of twosterically-stabilized colloid species. Such systems are of-ten modelled as a highly size-asymmetric binary mixtureof hard spheres, for which the effective interactions arisefrom the celebrated depletion mechanism . As well asbeing a key ingredient in determining effective Hamilto-nians for asymmetrical hard sphere mixtures, depletioninteractions can be directly measured in experiments .For a pair of big hard spheres in a sea of small hardspheres, the effective potential takes the form φ eff ( r bb ) = φ bb ( r bb ) + W ( r bb ) , (1)where φ bb ( r bb ) is the bare hard sphere potential be-tween two big spheres of diameter σ b whose centers areseparated by a distance r bb , and W is the “depletionpotential” which is mediated by the small spheres ofdiameter σ s . In this paper we consider additive hardsphere mixtures so that the big-small interaction diame-ter σ bs = ( σ b + σ s ) /
2. In that case the depletion poten-tial is attractive for small separations of the big spheres,but decays in an exponentially damped oscillatory fash- ion at large separations. The physics of the attraction iswell understood: The exclusion or depletion of the smallspheres as the big ones come close together results inan increase in free volume available to the small speciesleading to a net increase of entropy .A number of theoretical prescriptions exist for deter-mining effective potentials, including integral equations(as summarized in the recent article by Bot¸an et al ),DFT (see the summary in Ashton et al ) and morpho-metric theory. However, these theoretical treatmentsinvolve approximations, the validity of which need tobe checked. Computer simulation potentially provides aroute to estimating effective potentials which is in princi-ple exact, and can therefore be used to verify theoreticalpredictions. Unfortunately, it too finds the regime oflarge size asymmetry extremely challenging. The diffi-culty stems from the slow relaxation of the big particlescaused by the presence of the small ones. Specifically, inorder to relax, a big particle must diffuse a distance oforder its own diameter σ b . However, for small size ra-tios, q ≡ σ s /σ b and even at quite low volume fractions ofsmall particles, very many small particles will typicallyoccupy the space surrounding a big particle and thesehem it in, greatly hindering its movement. In compu-tational terms this issue mandates a very small Molecu-lar Dynamics timestep in order to control integration er-rors, while in basic Monte Carlo (MC), a very small trialstep-size must be used in order to maintain a reasonableacceptance rate. Consequently, the computational costof simulating highly size asymmetric mixtures by tradi-tional means is prohibitive at all but very low volumefractions of small particles.Owing to these difficulties, most previous simulationstudies of highly size asymmetrical ( q (cid:46) .
1) hard spheremixtures have adopted an indirect route to measuringdepletion potentials based on measurements of interpar-ticle force . The strategy rests on the observation that theforce between two big particles can be expressed in termsof the contact density of small particles at the surface ofthe big ones . By measuring this (angularly depen- a r X i v : . [ c ond - m a t . s o f t ] J un dent) contact density for fixed separation r bb of the bigparticles and repeating for separations ranging from con-tact, r bb = σ b , to r bb = ∞ , one obtains the force profile.This can in turn be integrated to yield an estimate ofthe depletion potential. However, the statistical qual-ity of the data obtained via this route is typically quitelow, particularly at small q ≤ . q ≤ . for a hard sphere study and Luijtenand coworkers for more general potentials. Thesestudies deployed a cluster algorithm (to be described insec. V A) to deal with the problem of slow relaxationoutlined above. However, this algorithm is limited in therange of particle volume fractions for which it will op-erate efficiently and thus there is a need for alternativeapproaches that extend this range to higher values.An additional drawback of previous studies is thatthey have treated the small particles canonically ratherthan grand canonically. Doing so complicates comparisonwith theoretical studies which are typically formulated interms of an infinite reservoir of small particles. It is alsoat variance with the common experimental situation of adepletant that is in equilibrium with a bulk reservoir.In what follows we consider how Monte Carlo simula-tion can be used to obtain direct and accurate estimatesof the depletion potential between two big hard spheresseparated by a distance r bb immersed in a dense sea ofsmall particles at size ratio q = 0 .
1. Our focus is onthe range of available techniques, their implementationand their relative utility; comparisons with theoreticalpredictions have appeared elsewhere . II. SYSTEM SETUP
The simulation setup that we consider for the measure-ment of depletion potentials is depicted in cross sectionin Fig. 1. It comprises a cuboidal periodic simulation boxwith dimensions L x = 3 . σ b , L y = L z = 2 . σ b . This boxaccommodates two big hard spheres and a large num-ber of small ones. (Though as described below, in someinstances it will prove beneficial to take one of the bigparticles to be a hard shell). Owing to the sphericalsymmetry of the depletion potential we can, without lossof generality, fix the center of one of the big particles atthe origin, while constraining the center of the other tooccupy points along the x -axis at x = r bb ≥ σ b . The onlyexception to this arrangement is the cluster algorithm tobe discussed separately in Sec. V A.We set the size of the small particles to be σ s = 0 . σ b ,ie. q = 0 .
1. We also elect to treat them grand canonicallyso that their total number fluctuates. Conceptually thiscorresponds to a colloidal system connected to a reser-
FIG. 1. A cross section through a snapshot of a configurationas described in the text. The simulation box contains a pairof big hard spheres, one of which is fixed at the origin, whilethe other is located at x = r bb , y = 0 , z = 0, with r bb = 1 . q = 0 .
1) at reservoir volumefraction η r s = 0 .
32. The section shown corresponds to theregion z < voir of depletant particles whose properties are param-eterized in terms of either the reservoir volume fraction η r s = πρ s σ /
6, (with ρ s = N s /V the reservoir numberdensity) or equivalently the conjugate chemical poten-tial µ rs . In practical terms, use of the grand canonicalensemble aids relaxation of small particle configurationsbecause particle transfers (insertions and deletions) canbe performed very efficiently. However to utilize this en-semble one needs to know accurately the chemical po-tential corresponding to a given η r s . We obtain thisfrom the equation of state of Kolafa et al , which wehave checked provides a highly accurate representation ofgrand canonical ensemble simulation data. Transfers ofsmall particles are effected using a standard grand canon-ical approach . For the most part we consider the caseof a rather high reservoir volume fraction of small parti-cles, η r s = 0 .
32, which also corresponds to the conditionsdepicted in the configurational snapshot of Fig. 1.
III. OVERVIEW OF COMPUTATIONAL STRATEGIES
We shall investigate two distinct routes to obtainingestimates of depletion potentials which we outline herebefore going into detail in Secs. IV and V. The first routeis based on measurements of the insertion probability ofone big sphere in the presence of the other; the second isbased on direct sampling of free energy differences asso-ciated with variations in the separation between the twobig spheres.
A. Insertion route and the shell trick
Let µ ex ( r bb ) be the excess chemical potential associ-ated with inserting a big sphere at some prescribed dis-tance r bb from another big sphere. It is straightforwardto show that this function is equivalent to the effectivepotential up to an additive constant i.e. W ( r bb ) = µ ex ( r bb ) − C, (2)where the constant C = lim r bb →∞ µ ex ( r bb ) . (3)To facilitate estimates of the excess chemical potential,one can appeal to the Widom insertion formula , whichin the case of hard particles reads µ ex ( r bb ) = − β − ln [ p i ( r bb )] . (4)Here p i ( r bb ) is the probability that an attempt to inserta big particle at x = r bb incurs no overlaps with smallparticles; it is calculated with respect to the ensemble ofconfigurations of the small particles. β is the inverse tem-perature, which in hard particle systems simply serves tobestow free energies with the appropriate dimensions; ac-cordingly we shall henceforth set it to unity.It follows from Eqs. 2-4 that the depletion potentialcan be expressed in terms of insertion probabilities as W ( r bb ) = ln (cid:18) p i ( ∞ ) p i ( r bb ) (cid:19) , (5)where p i ( ∞ ) represents the insertion probability for in-finite separation of the big spheres, which in practicalterms can be determined as the insertion probability ofa big sphere in a simulation box containing only smallparticles.The computational task is then to measure the inser-tion probability p i ( r bb ). Unfortunately, for the values of η r s of interest this probability is almost vanishingly small,a fact which renders simple sampling ineffective. Conse-quently we adopt a bespoke ‘gradual insertion’ approach,based on the use of tunable interactions and biased MonteCarlo sampling. Details of this approach are postponeduntil Sec. IV. Here it suffices to note that in implement-ing such schemes a very useful “geometrical shortcut”derives from the fact that it is not actually necessary toconsider the insertion probability of a big hard sphere in order to calculate the depletion potential. Instead itis sufficient and (generally much more efficient) to mea-sure the insertion probability for a hard shell of diameter σ b having infinitesimal thickness, as shown in the snap-shot of Fig. 2. The essential observation is that whenfully inserted, a hard shell particle encloses a number of small particles and although these remain in equilibriumwith the reservoir (by means of particle transfers) theyare fully screened from the rest of the system becausetheir surfaces cannot penetrate the shell wall. Thus thecontribution to the partition function from the enclosedparticles is independent of r bb , and therefore representsa constant contribution to µ ex ( r bb ) which vanishes fromthe difference in Eq. 2. Accordingly Eq. 5 applies equallyto shell insertion as it does to sphere insertion. Of coursefrom a computational standpoint, the task of inserting ahard shell is much less challenging than that of inserting ahard sphere (as can be appreciated by comparing Figs. 1and 2): essentially the insertion probability falls with theparticle size ratio like q rather than q . Shell insertion isdeployed in each of the three gradual insertion methodsto be described in Sec. IV. FIG. 2. A cross section through a configuration containing abig hard sphere fixed at the origin (left) and a fully insertedbig hard shell (right) in equilibrium with a fluid of small par-ticles at η r s = 0 .
32. The particle size ratio is q = 0 . r bb = 1 .
19. The section shown corresponds to z < A further geometrical shortcut results from noting thatthe convergence of the ensemble average over small parti-cle configurations required to calculate the shell insertionprobability depends on how quickly the small particles inthe region of the shell decorrelate. To enhance this re-laxation rate we preferentially perform grand canonicalinsertions and deletions of small particles within a shellsubvolume of radius 0 . σ b ≤ r ≤ . σ b centered on theshell. Updates inside the subvolume occur with a fre-quency 50-fold that of outside. This approach –whichsatisfies detailed balance– greatly reduces the time spentupdating small particles whose coordinates are relativelyunimportant for the quantity we wish to estimate. B. Direct sampling route
The particle insertion approach outlined above relieson extracting the depletion potential from differences inthe measured values of the insertion probability as a func-tion of r bb . However, even when using the shell insertiontrick, the difference ln p i ( ∞ ) − ln p i ( r bb ) that providesthe depletion potential via Eq. 5, is (notwithstanding thelogarithm) typically small compared to the absolute val-ues of ln p i ( ∞ ) and ln p i ( r bb ). Potentially, therefore, agreat deal of computational effort is required to obtain areasonable accuracy in W ( r bb ). In view of this we haveinvestigated an alternative strategy for obtaining the de-pletion potential which directly measures changes in thefree energy as the separation between the two big spheresis varied. To achieve this, however, specialist methodsare required to overcome the steric hindrance to the dis-placement of a big particle in a sea of much smaller ones.In section V we consider two methods that enable suchdisplacements via collective updates of a big sphere andmany small ones. They are: (i) the cluster algorithm ofDress and Krauth , which allows the depletion poten-tial to be built up directly from the sampled histogram ofbig particle separations, and (ii) a new constrained biasedcluster move, which permits estimates of the free energydifference associated with a prescribed displacement of abig particle. IV. INSERTION ROUTE: IMPLEMENTATIONS
In this section we outline three methods that exploitthe insertion route to determine the depletion potential.The basic idea is to to fix a hard sphere at the origin andthen estimate the probability of inserting a hard shell atcoordinates x = r bb , y = z = 0. In practice, however,for highly size asymmetrical mixtures and at all but thesmallest values of η r s , simple sampling of the insertionprobability is too inefficient to yield accurate results. In-stead a more elaborate gradual insertion technique is re-quired to render the approach feasible. We note thatkey elements of the relevant strategies and general sam-pling issues for determining insertion probabilities (andthence excess chemical potentials) have been discussedpreviously elsewhere , though not in the context ofhighly size asymmetrical fluid mixtures. A. Method I: Expanded ensemble
This method, which has been briefly reportedpreviously draws on earlier related studies. It in-volves defining an extended set of states for the inter-action between the shell particle and the small particlesand implementing Monte Carlo updates that make tran-sitions between these states.
1. Description
To estimate p i ( r bb ) for the shell we suppose that itcan exist in one of M possible ‘ghost’ states or ‘stages’in which it interacts with a small hard sphere (a distance r bs away) via the potential φ (m)g ( r bs ) = (cid:40) − ln λ (m) , ( σ b − σ s ) / < r bs < σ bs , otherwise. (6)Here m = 0 . . . M − ≤ λ (m) ≤ λ (m) > finite so that overlaps between small parti-cles and the big one can occur. If we denote by N o theinstantaneous number of such overlaps, then the config-urational energy associated with the shell in stage m isΦ (m)g = − N o ln λ (m) . (7)Clearly for λ (m) = 1, the shell is completely non-interacting, while for λ (m) = 0 it is infinitely repulsive.To span this range we set the extremal stages λ (0) = 1and λ ( M − = 0 (in fact we choose λ ( M − = 10 − toavoid numerical infinities), and define a set of M − λ (m) , m = 1 , . . . , M − , . . . , M −
1, i.e. that permits the shell interaction tofluctuate smoothly between the two extremes of interac-tion strength.Details of a suitable Metropolis scheme for samplingthe full range of m = 0 . . . M − The basic idea is to perform grand canon-ical simulation of the small particles, supplemented byMC updates that allow transitions m → m (cid:48) = m ± p a ( m → m (cid:48) ) = min (cid:16) , exp [ − (Φ (m (cid:48) )g − Φ (m)g ) + ∆ w ] (cid:17) , (8)where ∆ w = w (m (cid:48) ) − w (m) , with w (m) a prescribed weightassociated with stage m (see below). Note that for tran-sitions that depart from the extremal stages m = 0 orm = M −
1, it is necessary to reject proposals that wouldtake m outside the range (0 , M − M stagesrepeatedly, permitting a histogram ˜ H (m) of their rela-tive probabilities to be accumulated. From this biasedhistogram, one unfolds the weight factors to obtain anestimate of the unbiased histogram: H (m) = ˜ H (m) exp ( w (m) ) . (9)After normalizing to unit integrated weight, this his-togram provides an estimate of the relative probability p (m | r bb ) of finding the system in each of the M stages.The insertion probability is simply the relative probabil-ity of finding the system in the extremal stages: p i ( r bb ) = p ( M − | r bb ) p (0 | r bb ) , (10)from which the effective potential (up to a constant) fol-lows via Eq. 5. Repeating the measurement for a suc-cession of values of r bb allows construction of the entiredepletion potential.
2. Remarks and results
The implementation of method I entails a certain de-gree of preliminary work. Firstly one must decide on thenumber of stages M and their locations in λ ∈ [0 , { λ (m) } , m = 1 . . . M − λ (0) = 1 and λ ( M − = 10 − . It is important that these choices resultin MC transitions m → m ± λ ,and (in short runs) measure the distribution of overlaps p ( N o | λ (m) ) for each. From this set we select a subset of M stages for which the acceptance rate for transitionsm → m ± M , while a largeacceptance rate necessitates a correspondingly larger M .Although we find empirically that the overall efficiency ofthe method is not particularly sensitive to the choice ofacceptance rate (provided it lies in the range 10% − M stages and the transition rate.Secondly one needs to prescribe a suitable set of M weights { w (m) } for use in the acceptance probabilityEq. 8. The role of these weights is to bias the accep-tance rates such as to enhance the sampling of statesof low probability. Generally speaking a suitable set ofweights is one which ensures approximately uniform sam-pling of the M stages . The weights can be determinedusing a variety of methods, though we favor the Tran-sition Matrix Monte Carlo (TMMC) method detailed inAppendix A. Note that having determined a suitable setof weights for one value of the big particle separation r bb ,this set will (typically) perform adequately at all valuesof r bb to be studied, at least provided the variations inthe depletion potential are not too large, as is certainlythe case for the range η r s ≤ .
32 considered here. Sim-ilarly one does not have to choose a new set of { λ ( m ) } for each choice of the big particle separation r bb , a singlechoice performs adequately for all separations. Fig. 3 shows data accumulated for r bb = 1 . , q =0 . , η r s = 0 .
32. For this state point, M = 16 stages wererequired to realize a 20% acceptance rate for transitionsin m. A portion of the time series resulting from thesampling of m is shown in Fig. 3(a), giving an impres-sion of the timescale over which the sampling covers theentire range. The estimates of the probability distribu-tion p ( λ (m) | r bb ) that results from unfolding the weightsfrom the measured histogram ˜ H (m) (cf. Eqs. 9 and 10)is shown in Fig. 3(b). From this, the insertion probabil-ity can be read off directly; it is found to be O (10 − ),demonstrating the scale of the depths in probability thatthe method allows one to plumb. The rationale for theextreme improbability of successfully inserting a shellwithout the support of biased sampling is to be foundin Fig. 2, specifically in the tightness of the small parti-cle packing at this value of η r s .Finally in this subsection we remark that since the fulldepletion potential is built up from separate and indepen-dent measurements of the insertion probability at variousvalues of r bb , there is the opportunity to exploit paral-lelism by farming out each measurement on multi-coreprocessors. B. Method II: Multiple overlapping histograms
Our second approach is related to the previous one inthat a set of M stages are used to control the strengthof interaction between the shell and the small parti-cles in the manner described by Eq. 6. The differ-ence is that here we don’t actually implement transitions λ (m) → λ (m ± , instead we simply measure the free en-ergy difference between successive values of λ via an exactfree energy perturbation method.
1. Description
The relevant expression for calculating free energy dif-ferences is the well known formula of Zwanzig , whichin our case, for a transition m → m (cid:48) = m + 1 reads: F (m (cid:48) ) − F (m) = − ln (cid:68) exp (cid:104) − (Φ (m (cid:48) ) g − Φ (m) g ) (cid:105)(cid:69) m ,r bb = − ln (cid:28) exp (cid:20) N o ln λ (m) λ (m (cid:48) ) (cid:21)(cid:29) m ,r bb = − ln (cid:32)(cid:88) N o (cid:18) λ (m) λ (m (cid:48) ) (cid:19) N o P ( N o | λ (m) , r bb ) (cid:33) . (11)Here the ensemble average is with respect to the smallparticle configurations in stage m, given a big particleseparation r bb .We can apply this formula in the forward and reversedirections, averaging the result to find: CPU Time (seconds) m (a) m -300 -240 -180 -120 -60 p(m) r bb =1.00 r bb =1.06 r bb = ¥
14 15 16 m -200 -196 (b) FIG. 3. (a) λ (m) vs CPU time on a 2 GHz processor for r bb = 1 . η r s = 0 . , q = 0 . M = 16 stages, but constitutes onlya small portion of the full run which comprised 35 CPUhours. (b) the unfolded histogram p ( λ (m) | r bb ) at a selec-tion of values of r bb . Differences in the insertion probabil-ity p ( λ ( M − | r bb ) /p ( λ (0) | r bb ) (inset) provide estimates for thevariations in the depletion potential. F (m (cid:48) ) − F (m) = 12 ln (cid:80) N o (cid:16) λ (m (cid:48) ) λ (m) (cid:17) N o P ( N o | λ (m (cid:48) ) , r bb ) (cid:80) N o (cid:16) λ (m) λ (m (cid:48) ) (cid:17) N o P ( N o | λ (m) , r bb ) . (12)Thus, operationally, having chosen a suitable set of in-termediates { λ (m) } , one simply measures the distributionof overlaps P ( N o | λ (m) , r bb ) at each λ (m) . This yields theinsertion probability vialn p i ( r bb ) = F (0) − F ( M − . (13)The depletion potential then follows by repeating thismeasurement for a sequence of values of r bb and utilizing Eq. 5 as was done in Sec. IV A.
2. Remarks and results
For this method to yield accurate results, stages haveto be placed at appropriate values of λ such that succes-sive distributions p ( N o | λ i ) and p ( N o | λ i +1 ) overlap sig-nificantly. This is essentially the same criteria for choos-ing the set of intermediates { λ (m) } that is required toyield a reasonable acceptance rate between all stages inmethod I (Sec. IV A). Indeed comparing with Eq. 8, onesees that Eq. 11 provides a measure of the acceptancerate for transitions between neighbouring stages as ex-plicitly implemented in method I. Accordingly it servesas a basis for thinning out, appropriately, the trial setof 1000 stages as described in Sec. IV A 2. The resultingset { λ (m) } is then equally applicable to methods I andII. We emphasize that for either method there is no needto recalculate the set { λ (m) } for each r bb of interest; de-termining a set for one value of r bb suffices for all valuesprovided the depletion potential does not vary by morethan a few k B T . We also remark in passing that whilemethod II bears some resemblance to thermodynamic in-tegration schemes , the estimates of the free energy dif-ferences are in principle exact- no numerical quadratureis involved.Fig. 4 shows our measurements of the set of M = 16individual distributions p ( N o | λ (m) ) for η r s = 0 . , q = 0 . { λ (m) } is thesame as that used in method I and is listed in the key. N p(N ) m 0 1 2 3 4 5 6 7 ln( ) 0.0 0.36 0.728 1.104 1.464 1.832 2.208 2.584m 8 9 10 11 12 13 14 15 ln( ) 2.992 3.392 3.848 4.368 5.0 5.856 7.24 FIG. 4. The measured form of the overlapping distributions p ( N o | λ (m) ) employed to measure the shell insertion proba-bility for r bb = 1 . , η r s = 0 . , q = 0 . λ are shownin the key. The distribution for m = 15 is not depicted as itencompasses only the N o = 0 state. The chief merit of the multiple overlapping histogramapproach compared to the expanded ensemble approach(method I) is its simplicity: no weights need to be calcu-lated before one can start to accumulate data. Its maindisadvantage compared to method I, is the need to per-form M independent simulations and synthesize the re-sults in a pairwise fashion. However, this drawback issomewhat mitigated by the fact that the independenceof the simulations for each λ (m) renders them triviallyparallel. Accordingly, one can farm out the calculationsfor each to a separate processor on a multiprocessor com-puter. Similarly the estimates of the insertion probabilityat the various values of r bb that are needed to constructthe full depletion potential are also independent, and cantherefore be accumulated in parallel. C. Method III: Umbrella sampling
This approach, which has some commonality with theumbrella sampling approach of Ding and Valleau , isconceptually simpler than the previous two in that it dis-penses with staged intermediates.
1. Description
The algorithm considers an imaginary shell of diame-ter σ b centered on x = r bb . The instantaneous numberof small particles, N o , that overlap this notional shellfluctuates with time, and hence one can measure its dis-tribution p ( N o | r bb ) as a histogram. Typically N o will belarge, but we can performs biased (“Umbrella”) samplingwith respect to insertion and deletion of the small par-ticles in order to accurately measure the probability ofstates having N o = 0. A little thought shows that thisprobability is just the shell insertion probability requiredfor Eq. 5.Operationally, transfers of small particles are per-formed according to the biased acceptance probabilities: p a ( N s → N s +1) = min (cid:18) , VN s +1 e µ + W + (cid:19) ,p a ( N s → N s −
1) = min (cid:18) , N s V e − µ + W − (cid:19) . (14)These are the standard criteria for the grand canonicalensemble , modified by a weight factor W ± that is non-zero if the proposed insertion or deletion of a small par-ticle leads to a change in the number of overlaps N o .Specifically W + = w (cid:0) N o ( { r } N s +1 ) (cid:1) − w (cid:0) N o ( { r } N s ) (cid:1) ,W − = w (cid:0) N o ( { r } N s − ) (cid:1) − w (cid:0) N o ( { r } N s ) (cid:1) . (15) Here N o ( { r } N s ) is the number of overlap arising from theset of position vectors { r } N s = r , r . . . r N s of N s smallparticles, while w ( N o ) is a weight function defined on thenumber of overlaps. These weights allow a single simu-lation run to sample not just the values of N o that aretypical for a given η r s , but also the entire range down to N o = 0. Accordingly one can measure a histogram of theweighted probabilities ˜ H ( N o | r bb ), from which the Boltz-mann histogram is obtained by unfolding the weights: H ( N o | r bb ) = ˜ H ( N o | r bb ) e w ( N o ) . (16)After normalization, this yields the probability distribu-tion p ( N o | r bb ), from which the insertion probability isread off as p ( N o = 0 | r bb ). The depletion potential (up toa constant) follows via eq. 5. Repeating for a sequence ofvalues of r bb allows one to build up the entire depletionpotential.
2. Remarks and results
As with method I, an appropriate set of weightsis required for this method to operate effectively andagain these can be readily determined using the TMMCmethod (Appendix A). Fig. 5 shows a time series ofthe sampled values of N o that results once the weightsare in place. Owing to the biasing, the system samplessmoothly the entire range from the most probable num-ber of overlaps N o ≈ N o = 0. Theresulting form for p ( N o | r bb ), obtained by unfolding theeffects of the weights and normalizing the resulting his-togram is shown in Fig. 5(b). From this one simply readsoff the shell insertion probability as p ( N o = 0 | r bb ).The chief merit of method III compared to methods Iand II is that it is parameter free: there are no stagedintermediates and therefore the associated inconvenienceand startup costs of determining their number and ap-propriate placement are obviated. Nevertheless the com-putational cost of calculating weights represents a signif-icant overhead as will be discussed in Sec. VI. We notethat method III is parallelisable, but only with respectto the separate measurements at various r bb needed tobuild up the depletion potential. V. DIRECT SAMPLING ROUTE: IMPLEMENTATIONS
We now turn to consider two schemes that accumulatethe depletion potential by focusing on the difference in ef-fective potential as one varies r bb . They both rely on col-lective (cluster) updates of big and small particles. Oneis based on the cluster algorithm of Dress and Krauth ,the other is a bespoke constrained cluster algorithm. CPU time (seconds) N (a) N o -250 -200 -150 -100 -50 p(N o ) r bb =1.00r bb =1.06r bb = ¥ o -205 -200 -195 P ( N o ) r bb =1.0r bb =1.06r bb = ¥ (b) FIG. 5. (a) N o ( t ) vs CPU time on a 2 GHz processor at r bb =1 . σ b , η s = 0 . , q = 0 . N o = ¯ N o to N o = 0, but constitutesonly a small portion of the full run which comprised 35 CPUhours. (b) The form of the overlap probability distribution p ( N o | r bb ) at three values of r bb . Differences in p ( N o = 0 | r bb )as a function of r bb (inset) yield the depletion potential asdescribed in the text. A. Method IV: Geometrical Cluster Algorithm
An efficient cluster algorithm capable of dealing withhard spheres mixtures was introduced by Dress andKrauth in 1995 . It was subsequently generalized toarbitrary interaction potentials by Liu and Luijten who dubbed their method the Geometrical Cluster Algo-rithm (GCA). A restricted Gibbs ensemble version of theGCA suitable for studying phase transitions was also sub-sequently developed . Here we describe the GCA fora general system of hard spheres in the canonical ensem-ble, before specializing to the case of a size asymmetricalbinary mixture.
1. Description
The particles comprising the system are assumed tobe contained in a periodically replicated cubic simula-tion box of volume V . The configuration space of theseparticles is explored via cluster updates, in which a sub-set of the particles known as the “cluster” is displacedvia a point reflection operation in a randomly chosenpivot point. The cluster generally comprises both bigand small particles and by virtue of the symmetry ofthe point reflection, members of the cluster retain theirrelative positions under the cluster move. Importantly,cluster moves are rejection-free even for arbitrary inter-particle interactions . This is because the manner inwhich a cluster is built ensures that the new configura-tion is automatically Boltzmann distributed.For hard spheres (there is no advantage in using shellsin this context), the cluster is constructed as follows: oneof the particles is chosen at random to be the seed parti-cle of the cluster. This particle is point-reflected with re-spect to the pivot from its original position to a new posi-tion. However, in its new position, the seed particle mayoverlap with other particles. The identities of all suchoverlapping particles are recorded in a list or “stack”.One then takes the top-most particle off the stack, andreflects its position with respect to the pivot. Any par-ticles which overlap with this particle at its destinationsite are then added to the bottom of the stack. This pro-cess is repeated iteratively until the stack is empty andthere are no more overlaps.Note that cluster updates only displace particles, theydo not allow their number to fluctuate. Accordingly, inorder to treat the small particles grand canonically, wealso perform insertions and deletions of small particleswith a chemical potential corresponding to the prescribed η r s , as outlined in Sec. II,The effective potential W ( r ) between two big parti-cles is defined in terms of the radial distribution function g ( r bb ), measured in the limit of infinite dilution: W ( r ) = − lim ρ b → ln[ g ( r bb )] , (17)for r bb > σ b . In our simulation studies this limit is ap-proximated by placing a single pair of big hard spheresin the simulation box. A finite-size estimate to g ( r bb ),which we shall denote g L ( r bb ), is then obtained by fixingthe first of these particles at the origin and measuring (inthe form of a histogram) the probability p ( r bb ) of findingthe second big particle in a shell of radius r bb → r bb + dr .Then g L ( r bb ) = p ( r bb ) p ig ( r bb ) , (18)where the normalization relates to the probability of find-ing an ideal gas particle at this radius: p ig ( r bb ) = 4 πr V . (19)To effect the measurement of g L ( r bb ), we modify theGCA slightly as follows: we choose one big particle to bethe seed particle, which we place randomly within a shell σ b < r bb < L/
2, centered on the second big particle,with L the linear box dimension. The location of thepivot is then inferred from the old and new positions ofthe seed particle. Thereafter clusters are built in thestandard way. This strategy ensures that we efficientlysample separations of the big particles that lie in therange σ b < r bb < L/ g ( r bb ) can sensibly bedefined for hard spheres in a cubic box.
2. Remarks and results
For the systems of interest in this work, we find thatthe GCA is efficient for reservoir packing fractions η r s ≤ .
2. Above this value, practically all the particles join thecluster, which merely results in a trivial point reflectionof the entire system. Indeed the efficiency drop is soprecipitous that η r s = 0 . . Doing so has been reported to extend theoperating limit to η r s (cid:39) .
34. However, for the case ofhighly asymmetrical mixtures we find that this strategydoes not significantly decrease the number of particlesin the cluster because as soon as a big particle joins thecluster and is point reflected it causes many overlaps withsmall particles.Fig. 6(a) shows the measured form of g L ( r bb ) for η r s =0 . , q = 0 . V = (3 σ b ) . For this measurement to providean estimate of W ( r bb ), it first has to be corrected forfinite-size effects, manifest in the failure of the functionto approach unity at large r bb . This is done (as has alsobeen described elsewhere ) by measuring the cumulativeintegral G ( R ) = (cid:90) R g L ( r ) dr . (20)This integral tends towards a smooth linear form quiterapidly as the upper limit R increases. The measuredlimiting gradient, ξ , of G ( R ) provides the requisite cor-rection factor according to g ( r ) = ξ − g L ( r bb ). FollowingEq. 17, the negative of the logarithm of g ( r ) then yieldsan estimate for the effective potential W ( r bb ), which isshown in Fig. 6(b).The most attractive feature of the GCA for determin-ing depletion potentials is that it allows direct sampling r g(r) r -4-3-2-101 W (r) FIG. 6. The measured form of g L ( r bb ) corresponding to η r s =0 . , q = 0 .
1, obtained via method IV for a simulation boxof dimensions V = (3 σ b ) . The limiting value differs fromunity due to the finite-size effects described in Ashton et al .The inset shows the depletion potential W ( r ) obtained byimplementing the finite-size correction described in the textto g L ( r bb ) and applying Eq. 17. of the quantity of interest without the need for multi-ple simulations or biased sampling. Its principal draw-back is that the method becomes unusable for η r s (cid:38) . η r s . The method described in the fol-lowing subsection achieves this, albeit at the expense ofintroducing biased sampling. B. Method V: Constrained cluster algorithm
In common with the GCA, this method collectivelymoves a big hard sphere and a number of small ones via aself inverse operation. However, in contrast to the GCA itis a constrained scheme in the sense that it measures thefree energy differences between two neighbouring discretevalues of r bb .
1. Description
The operation of the method is shown schematically inFig. 7. One big hard sphere (particle A ) is fixed at theorigin. The other (particle B ) can occupy discrete valuesof r bb set out along a one-dimensional radial grid whichwe take to be the x -axis. Let us label the grid pointsby the index i , and consider the situation when the bigparticle is stationed at r bb = x i . We then estimate thefree energy difference between grid points i and i + 1 inthe following manner.With particle B stationed at grid point i , we equili-brate the small particles via transfers with the reservoir.0 reflection plane N o =3 overlaps (a)(b) BAA B
FIG. 7. Schematic illustration of the constrained cluster up-date of method V. (a)
The big particle B at x i undergoesa plane reflection to x i +1 and thereby overlaps with n smallparticles which are themselves subsequently reflected in theplane. (b) In their new position, the n small particles overlapswith N o other particles (one of which may be the big particle A at the origin). The number of such secondary overlaps N o is the primary observable. For some equilibrium configuration of the small particleswe then consider (but do not implement) a trial move totake particle B from grid point i to grid point i + 1 asfollows:1. Reflect the center of particle B in the plane normalto the x axis which cuts the x axis at x = ( x i + x i +1 ) /
2. This takes particle B from grid point i togrid point i + 1 as shown in fig. 7(a).2. Under this move, particle B will overlap with anumber n , say, of small particles. We then imaginereflecting these n small particles in the same reflec-tion plane. This switches them into the space leftby particle B , see fig. 7(a).3. After undergoing this reflection, some of the n smallparticles will overlap with other small particles orwith the big particle A , as shown in fig. 7(b). Thenumber of such ‘secondary’ overlaps is the observ-able N o for the current configuration of small par-ticles.One then samples the fluctuations in N o with respectto the ensemble of small particle configurations and accu-mulates its probability distribution p ( N o ) as a histogram. Similarly to methods I-III, it is beneficial to preferen-tially implement transfers of small particles in a shellregion around big particle B ; this concentrates the com-putational effort on those regions which contribute mostto the measurement. The sampling of the small parti-cle configurations is biased so as to enhance the occur-rence of values of N o down to N o =0. This is achieved bydefining a weight function w ( N o ) which is incorporatedin the GCE acceptance probabilities Eq. 14, in exactlythe same manner as described for method III. An appro-priate weight function can be found automatically usingthe TMMC method described in the appendix.Together these measures enable an efficient and accu-rate estimate for the probability that the trial collectivemove leads to N o = 0, ie. a valid hard sphere configu-ration. Let us denote this probability p + i (0) because wehave measured it with particle B moving from grid point i to i +1. Similarly we can measure the probability p − i +1 (0)that a move from i + 1 → i leads to zero overlaps. Thenthe measured ratio p + i (0) /p − i +1 (0) provides the differencein the depletion potential between grid points i and i + 1via an an expression akin to Bennett’s acceptance ratioformula : W ( x i +1 ) − W ( x i ) = ln p − i +1 (0) p + i (0) . (21)From measurements of the difference in the depletionpotential between all neighbouring pairs of grid points,one extracts the depletion potential itself simply by sum-ming, commencing at a value of r bb sufficiently large that W ( r bb ) can be considered to have decayed to zero.
2. Remarks and results
Compared to the GCA (method IV), the principal as-set of method V is that it permits study of considerablylarger volume fractions of the small particles. This is be-cause the number of particles involved in the collectivemove is not allowed to grow indefinitely. Instead clus-ter growth is truncated after one iteration and biasedsampling used to obtain the information required to esti-mate the depletion potential. We note that a constrainedcluster algorithm suitable for estimating depletion po-tentials has previously been described by Malherbe andKrauth , however it does not truncate cluster growthand therefore is limited to much lower values of η r s thanthe present approach.In common with the gradual insertion methods I andIII, the constrained cluster algorithm requires (in gen-eral) knowledge of a set of weights for its operation. How-ever, because the method focuses on free energy differ-ences, the typical number of overlaps N o is generally farfewer than encountered in methods I and III, and hencethe degree of weighting required to reach N o = 0 is muchless. For example, for η r s = 0 .
32, and a grid point sepa-ration of x i +1 − x i = 0 .
05 we find N o ≈
20 (see Fig. 8)1which is to be compared with the ≈
200 overlaps that oc-cur for shell insertion in methods I and III. Thus weightcalculation is relatively quick and easy for method V, andindeed we find that if we reduce the small particle vol-ume fraction to η r s (cid:46) .
2, then no weights are required atall since the system samples the N o = 0 state sufficientlyoften without the aid of biasing. Even in cases whereweighting is required, it is in general not necessary tocalculate weights for every grid point; to the extent thatthe effective potential does not vary strongly betweengrid points, weights found for one grid point will sufficefor all other grid points. N o -6 -5 -4 -3 -2 -1 p i +1 ( N o ) p + i ( N o ) FIG. 8. Estimates of p + i ( N o ) and p − i +1 ( N o ) for i = 0, cor-responding to contact of the big particles (i.e. r bb = 1 . σ b )as obtained using method V. The grid point separation is x i +1 − x i = 0 .
05. The ratio of the values of these functionsfor N o = 0 provides an estimate of the difference in the effec-tive potential between the grid points via Eq. 21. Although it is perhaps reminiscent of methods that ob-tain the depletion potential by integrating the measuredforce in an MD setting , method V provides exact dif-ferences in the depletion potential i.e. no quadrature isrequired. However, one downside of the need to sumfree energy differences to obtain the depletion potentialis that cumulative errors arise. The error grows with thenumber of differences summed and can potentially leadto an estimate for W ( r bb ) that whilst appearing quitesmooth, nevertheless deviates significantly from the ex-act form. Since we commence summing the free energydifferences at large r , where the potential can be assumedto be essentially zero, this implies that the largest errorsoccur near contact. To be more precise, for j = 1 · · · N free energy differences, the variance in the sum is simplythe sum of the variances of the individual (uncorrelated)estimates ie σ N = (cid:80) Nj =1 σ j . If each individual measure-ments receives an equal computational expenditure thento a good approximation the cumulative error after sum-ming N differences is simply σ N = √ N σ . This growth inthe uncertainty in the estimate of W ( r ) as r decreases,contrasts with the gradual insertion methods where everypoint in the estimate of W ( r bb ) is independent. Finally in this subsection we remark that in commonwith methods I-III, method V is parallelisable with re-spect to calculations along the grid: one can simply setup independent copies of the simulated system each ofwhich calculates W ( x i +1 ) − W ( x i ) for a different gridpoint i . VI. DISCUSSION
In the preceding two sections we have described fivedistinct methods for determining depletion potentials inhighly size asymmetrical hard sphere mixtures. We nowturn to a discussion of their relative merits.Let is begin by comparing the gradual insertion meth-ods I-III amongst themselves. In terms of their relativeefficiency, we find that once prepared so that samplingcan commence, each of the methods I-III take a simi-lar amount of CPU time to achieve a given statisticalaccuracy for W ( r bb ). This is shown in Fig. 9(a) whichdisplays the form of the depletion potential at η s = 0 . q = 0 . r bb , renders it slightly inferior to method III, in ourview. That said, and at the end of the day, whether onechooses to use method II or III is probably as much amatter of personal taste than of efficiency.In terms of the domain of applicability of the gradualinsertion approach, methods I-III, we note that all threemethods are effective in facilitating estimates of depletionpotentials at rather high volume fractions of small parti-cles. In this paper we have presented results for systemshaving η r s (cid:46) .
32 and the size ratio q = 0 .
1. Elsewhere we have shown that gradual shell insertion operates ef-2 r/ σ b -6-4-202 W ( r ) Method IMethod IIMethod III r/ σ b -6-4-2024 W ( r ) Method IIIMethod V
FIG. 9. Estimates of the form of W ( r bb ) for q = 0 . , η rs =0 . (a) Comparison of the results of methods I-III. (b)
Comparison of the results of methods I and III. All data cor-responds to a computational investment of 35 CPU hours perpoints on a 2 GHz processor. For methods I-III statisticalerrors are comparable with symbol sizes. For method V, thestatistical error accumulates in going from large to small r asindicated by the representative error bars. fectively up to about η r s = 0 .
35. This limit arises from arapid increase in the relaxation time for the small parti-cles which are tightly packed at this volume fraction. Forsize ratios smaller than q = 0 .
1, the problem of determin-ing the depletion potential is certainly computationallyharder than for q = 0 . N o between the shell and the small particles isgreater. Nonetheless we still expect, in principle, to beable to reach small particle volume fractions of η r s ≈ . W ( r bb ) across a set of values of r bb36 .Since each such measurement is independent, this has theattractive feature that there are no correlations amongthe data points that form the estimate of W ( r bb ). How-ever, a potential disadvantage of the approach arises fromthe fact that W ( r bb ) is obtained as the difference of twomeasurements, ie. W ( r bb ) = ln p i ( ∞ ) − ln p i ( r bb ). Ingeneral, both ln p i ( ∞ ) and ln p i ( r bb ) are large comparedto their difference, and thus, in effect, the gradual in- sertion approach calculates a small number by subtract-ing measurements of two large ones. Accordingly for agiven fractional uncertainty in ln p i ( r bb ), the correspond-ing fractional uncertainty in W ( r bb ) is larger by a fac-tor of √ p i ( r bb ) /W ( r bb ), requiring a greater compu-tational effort to obtain a satisfactorily smooth estimateof W ( r bb ). To give this issue some scale, for the case q = 0 . , η r s = 0 . − ln p i ( r bb ) ≈
200 for shellinsertion (while incidently, − ln p i ( r bb ) ≈
600 for sphereinsertion). These are to be compared with the maximumvariation in W ( r bb ) of <
4. Increasing η r s to 0 .
32 givesfor shell and sphere insertion, values of − ln p i ( r bb ) ≈ − ln p i ( r bb ) ≈ W ( r bb ) of ≈ η r s (cid:46) . η r s as large as those accessible to the gradual in-sertion methods. However, a caveat is that the apparentsmoothness of the estimates of W ( r ) arising from methodV may belie the true absolute error in W ( r ), which accu-mulates from large to small values of r . Our tests in thehigh density regime (cf. Fig. 9(b)), show that for a givenexpenditure of computational effort the maximum sta-tistical error in the potential obtained from method V iscomparable, but not significantly superior to the gradualinsertions methods. However, there is scope for furtherimproving the efficiency of method V by using a pair ofspherical caps rather than a spherical shell for the sub-volume in which preferential updating of small particlesis performed. Such subvolumes would include a higherproportion of the small particles that are effected by thevirtual move and thus increase the rate of fluctuation in N o VII. SUMMARY AND OUTLOOK
In summary, we have investigated a number of simula-tion techniques that facilitate accurate measurements ofdepletion potentials in highly size asymmetrical mixturesof hard spheres. Two categories of approach were consid-ered: (i) gradual insertion and (ii) cluster methods. Inthe first category, three flavors of methods were describedall of which obtain the depletion potential via measure-ments of the insertion probability of a big sphere or shellin the neighbourhood of another big sphere. Once pre-pared so that sampling could begin, all three insertionmethods showed comparable efficiency. However, differ-3ence were found in the startup costs associated with fac-tors such as whether the respective methods require pre-calculation of staged intermediates and/or weight factors.The gradual insertion methods allows one to obtain de-pletion potentials for small particle volume fractions ofup to about η r s = 0 .
35. However, to reach this limit itis essential to employ the ‘geometrical shortcuts’ that wehave described, namely shell insertion and preferentialsampling of small particles in the neighbourhood of a bigone. We remark that gradual insertion techniques haverecently been extended to systems containing many bigparticles in a full grand canonical ensemble simulationscheme for highly size asymmetrical fluid mixtures .In the second category, two cluster algorithms wereconsidered: the Geometrical Cluster Algorithm and abespoke constrained cluster method. The GCA is veryefficient provided η r s ≤ .
2. The constrained cluster algo-rithm considerably extends the range of η rs for which de-pletion potentials can be calculated to at least η r s ≈ . . Forthe cluster methods, a version of the GCA suitable forarbitrary potentials is well known . The constrainedcluster method could similarly be easily extended to ar-bitrary interactions by considering the energy associatedwith the trial move and measuring the ratio of acceptanceprobabilities for the forward and reverse move. ACKNOWLEDGMENTS
This work was supported by EPSRC grantsEP/I036192 and GR/F047800 and the Visiting Post-graduate Scholar Programme of the University of Bath.Some of the simulations were performed on a computerfunded by the HEFCE infrastructure fund. VSGgratefully acknowledges the support of CSIC and a JAEprogram PhD fellowship from the Direcci´on Generalde Investigaci´on Cient´ıfica y T´ecnica under GrantNo. FIS2010-15502 and from the Direcci´on Generalde Universidades e Investigaci´on de la Comunidad deMadrid under Grant No. S2009/ESP/1691 and ProgramMODELICO-CM. We thank Rob Jack and Bob Evansfor useful conversations. C. N. Likos, Phys. Rep. , 267 (2001). L. Belloni, J. Phys. Condens. Matter , R549 (2000). H. N. W. Lekkerkerker and R. Tuinier,
Colloids and the Deple-tion Interactions , Vol. 833 of
Lecture Notes in Physics (Springer,Berlin / Heidelberg, 2011). J. C. Crocker, J. A. Matteo, A. D. Dinsmore, and A. G. Yodh,Phys. Rev. Lett. , 4352 (1999). V. Bot¸an, F. Pesth, T. Schilling, and M. Oettel, Phys. Rev. E , 061402 (2009). D.J. Ashton, N.B. Wilding, R. Roth, and R. Evans, Phys. Rev.E , 061136 (2011). M. Oettel, H. Hansen-Goos, P. Bryk, and R. Roth, Euro. Phys.Lett. , 36003 (2009). T. Biben, P. Bladon, and D. Frenkel, J. Phys: Condens. Matter , 10799 (1996). R. Dickman, P. Attard, and V. Simonian, J. Chem. Phys. ,205 (1997). B. G¨otzelmann et al. , Europhys. Lett. , 398 (1999). A. R. Herring and J. R. Henderson, Phys. Rev. E , 011402(2007). P. Attard, J. Chem. Phys. , 3083 (1989). J.G. Malherbe and S. Amokrane, Mol. Phys. , 355 (2001). J. Liu and E. Luijten, Phys. Rev. E , 066701 (2005). S. A. Barr and E. Luijten, Langmuir , 7152 (2006). J. Kolafa, S. Labik, and A. Malijevsky, Phys. Chem. Chem. Phys. , 2335 (2004). D. Frenkel and B. Smit,
Understanding Molecular Simulation (Academic, San Diego, 2002). B. M. Mladek and D. Frenkel, Soft Matter , 1450 (2011). B. Widom, J. Chem. Phys. , 2808 (1963). C Dress and W Krauth, J. Phys. A , L597 (1995). I. Nezbeda and J. Kolafa, Mol. Sim. , 391 (1991). P. Attard, J. Chem. Phys. , 2225 (1993). N. B. Wilding and M. Muller, J. Chem. Phys. , 4324 (1994). D.A. Kofke and P.T. Cummungs, Mol. Phys. , 973 (1997). A. D. Bruce and N. B. Wilding, Adv. Chem. Phys , 1 (2003). A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov, and P. N.Vorontsov-Velyaminov, J. Chem. Phys. , 1776 (1992). D. J. Ashton and N. B. Wilding, Mol. Phys. , 999 (2011). R. W. Zwanzig, J. Chem. Phys. , 1420 (1954). K. Ding and J.P. Valleau, J. Chem. Phys. , 3306 (1993). J. Liu and E. Luijten, Phys. Rev. Lett. , 035504 (2004). J. Liu, N.B. Wilding, and E. Luijten, Phys. Rev. Lett. , 115705(2006). D.J. Ashton, N.B. Wilding, and P. Sollich, J. Chem. Phys. ,074111 (2010). D.J. Ashton, J. Liu, E. Luijten, and N.B. Wilding, J. Chem.Phys. , 194102 (2010). C.H. Bennett, J. Comput. Phys. , 245 (1976). J. G. Malherbe and W. Krauth, Mol. Phys. , 2393 (2007). Note also that methods I-III permit direct estimates of the con-tact value of the depletion potential. This contrasts with meth-ods that obtain the depletion potential by integrating the force,which rely on extrapolation to estimate the contact value. F.G. Wang and D.P. Landau, Phys. Rev. E. , 056101 (2001). P. Virnau and M¨uller, J. Chem. Phys. , 10925 (2004). G.R. Smith and A.D. Bruce, J. Phys. A , 6623 (1995). J.R. Errington, J. Chem. Phys. , 3130 (2004). G.C. McNeil-Watson and N.B. Wilding, J. Chem. Phys. ,064504 (2006).
Appendix A: Transition Matrix Monte Carlo