Multilevel simulation of hard-sphere mixtures Open Access

A key challenge in the molecular simulation of soft matter is posed by the separation of length-scales between its microscopic description and the existence or emergence of mesoscopic structure. In such cases, one often relies on coarse-grained (CG) descriptions of the system that (approximately) integrate out microscopic degrees of freedom^1–3 to yield a tractable simplified model. Examples in the soft-matter context include polymers,^4,5 biomolecules,^6–8 and colloidal systems.^9,10 Such CG descriptions are essential for multi-scale modeling approaches.^11,12 However, they are not usually exact, and the associated coarse-graining errors are often difficult to assess.

Such CG models have been studied extensively in colloidal systems with depletion interactions.^9,13,14 The typical example is a mixture of relatively large colloidal particles with much smaller non-adsorbing polymers that generate effective attractions between the colloids. This can drive de-mixing, crystallization, or gelation, depending on the context. Model systems in this context include the Asakura–Oosawa (AO) model⁹ where the CG model can even be exact if the disparity between colloid and polymer radii is large enough. The theoretical tractability of the AO model arises from a simplified modeling assumption that polymer particles act as spheres that can interpenetrate.

Alternatively, one may consider a mixture where both the colloids and the depletant are modeled as hard spheres. From a theoretical perspective, this is an interesting model in its own right as it undergoes a fluid–fluid de-mixing phase separation for large enough size-disparities and concentrations.^10,15–17 This happens despite the lack of attractive forces between the particles in the model and can be attributed to geometric packing effects of the big and small spheres.

Direct simulation of such mixtures is very challenging because of the large number of small particles. Accurate CG models are also available in this context,¹⁰ but the CG representations are not exact: their errors can be detected by accurate computer simulation of the full (FG) mixtures. Hence, such models are natural testing grounds for theories and simulation methods associated with coarse-graining.

In this context, we recently developed a method^17,18 that links a CG description with the underlying fine-grained (FG) description. We call this the two-level method because the CG and FG models describe the same system with different levels of details. The method was validated by computations on the AO system,¹⁸ where it provided numerically exact results for the FG model, even in the regime where the CG description is not quantitatively accurate. The methodology was also applied to the hard sphere mixture,¹⁷ where it provided a quantitative analysis of the critical point associated with de-mixing of the large and small particles.

These previous results rely on the idea that properties of the FG model can be estimated in terms of some CG quantity, with an additive correction that accounts for the coarse-graining error. This is an importance sampling method, familiar in equilibrium statistical mechanics from free-energy perturbation theory,¹⁹ which involves reweighting between two thermodynamic ensembles. In the present context, the reweighting factors depend on the free energy of the small spheres, computed in a system where the large particles are held fixed. This free energy can be estimated by an annealing process based on Jarzynski’s equality^20–22 that slowly introduces small particles to fixed CG configurations.

In this paper, we present an extension of the two-level method that incorporates additional intermediate levels to improve the overall performance. Specifically, we introduce a step in the annealing process where small particles are partially inserted in regions close to big particles. Before finishing the small-particle insertion, we then replace weighted sets of configurations with unweighted ones, duplicating configurations with large weight and deleting ones with low weight. This resampling step allows us to make optimal use of the information available at the intermediate stage, focusing our subsequent computations on configurations that matter.

This general approach fits in the framework of sequential Monte Carlo (SMC).^23–25 Such algorithmic ideas have been successfully applied in applications across disciplines under various names, including population Monte Carlo²⁶ or the go-with-the-winners strategy.²⁷ Examples in computational physics include the pruned-enriched Rosenbluth method for polymers,²⁸ the cloning method for rare events,²⁹ and diffusion quantum Monte Carlo.³⁰ We combine the SMC method with an additional variance reduction strategy. Instead of estimating the FG average directly, we combine a CG estimate with estimates of subsequent level differences using the previous levels as control variate. This is the idea behind multilevel Monte Carlo methods.^31–33 The combination of a difference estimate with SMC has been previously investigated, for example, in Refs. 34–36. As in Ref. 17, we develop a general method alongside its application to highly size-asymmetric binary hard-sphere mixtures, which provide a challenging but well understood example to benchmark our algorithm.

This paper is organized as follows: In Sec. II, we introduce the hard-sphere mixture model. In Sec. III, we summarize the setup of the two-level method before presenting our extension to three (or more) levels. The three-level method requires an intermediate level for the hard-sphere mixture, whose details we discuss in Sec. IV. In Sec. V, we present a numerical test of the method and compare its performance against the two-level method, and in Sec. VI, we present convergence results. We conclude in Sec. VII.

Throughout this work, we illustrate the multilevel method with an example system, which is a mixture of large and small hard spheres at size ratio 10:1. This system is challenging for simulation because the big particles may display interesting collective behavior (in particular, a critical point), but the dominant computational cost for simulation comes from the large number of small particles.

However, despite our focus on this single example, we emphasize that the multilevel method is presented in a general way, which should also be applicable in other systems with a separation of length scales.

The example system is a mixture of big and small particles, whose diameters are σ_B and σ_S, respectively. We consider a periodic box [0, L]³ of linear size L and we work in the grand canonical ensemble. (This choice is particularly relevant for analysis of de-mixing, where the number density of large particles is a suitable order parameter.³⁷)

In a given configuration, the numbers of big and small particles are N and n, respectively; the position of the ith big particles is R_i, while the position of the jth small particle is r_j. We denote the configurations of big and small particles by $C=(N;R1,…,RN)$ and $F=(n;r1,…,rn)$⁠, respectively, and the full configuration is denoted $X=(C,F)$⁠.

Since the particles are hard, the temperature plays no role in the following, and so we set the temperature as k_BT = 1 without any loss of generality. The equilibrium distribution of the mixture is described by a probability density

where the subscript F indicates that we refer to the FG model, μ_B and μ_S are the chemical potentials for the large and small particles, respectively, and Ξ_F is the grand canonical partition function. The particles are hard (non-overlapping), so the potential energy is

This p_F is normalized as $1=∫pF(C,F)dCdF$⁠, the precise meaning of these integrals is given in Appendix A 1.

Within this setting, the dimensionless parameters of the model are the ratio of the particle diameters σ_B/σ_S, the system size parameter L/σ_B, and the two chemical potentials μ_S and μ_B. In practice, μ_S is more naturally parameterized by the associated reservoir volume fraction $ηSr$⁠, which we relate to μ_B via an accurate equation of state.³⁸ 

Our multi-level method is designed for accurate estimates of properties of the large particles. Specifically, we consider observable quantities of interest $A=A(C)$ that only depend on the large particles. (Examples are discussed in Sec. II C: see also Fig. 3.) Our aim is to compute the equilibrium average of A, that is,

Since $A(C)$ does not depend on $F$⁠, it is natural to define the marginal distribution for the big particles

so that $⟨A⟩F=∫A(C)pF(C)dC$⁠. A similar situation occurs in the context of statistics, where one seeks to analyze the behavior of a few quantities of interest in a high-dimensional system: in that context, the small-particle degrees of freedom in (3) would be referred to as nuisance parameters. This means that their values are not required to compute the quantity of interest, but their statistical properties strongly affect the average of this quantity.

If samples for the marginal distribution $pF(C)$ could be generated by an MC method for the big particles alone, this would make the system much more tractable by simulation. This is a central idea in coarse-grained modeling.² However, the complexity of packing of the small hard spheres means that $pF(C)$ is a complex distribution, and it is not possible to sample it exactly. A great deal of effort has gone into developing CG models that approximate this distribution with high accuracy.^10,16,39

A suitable CG model is an equilibrium distribution with probability density

where Ξ_C is the partition function, and the CG (effective) interaction energy is

where V₂ is a pairwise interaction potential. Averages with respect to the CG model are denoted as

For a suitably chosen V₂, the coarse distribution $pC(C)$ can be an accurate approximation to $pF(C)$⁠. For the CG model in this work, we take the accurate potential V₂ = V_RED, developed by Roth, Evans, and Dietrich.¹⁰ Following Ref. 17, we choose Δμ such that the distributions of N coincide for FG and CG models.

Throughout the paper, we benchmark our numerical methods by considering the hard-sphere mixture with fixed parameters, as follows: We take the ratio of particle sizes (σ_B/σ_S) = 10, the linear size of the periodic system is L = 31σ_S, and the small-particle (reservoir) volume fraction is $ηSr=0.2$⁠. This volume fraction is large enough to generate a significant depletion attraction between the large particles, but not strong enough to cause de-mixing of the large and small particles.¹⁷ 

Aspects of the CG and FG models are illustrated in Figs. 1(a)–1(c) for these parameters. In particular, we show representative configurations of the CG and FG models, as well as a plot of the RED potential. While direct grand canonical Monte Carlo (GCMC) sampling of the full mixture is possible, in principle, it should be apparent from Fig. 1(c) that this would be intractable because insertion of large particles in such a fluid is hardly ever possible. Advanced MC methods^39,40 might be applicable, but these tend to struggle when the volume fraction gets large. This motivates the development of two-level and multi-level methods.

Figure 2 highlights properties of the distribution of the number of big particles for the CG model when varying the effective large-particle chemical potential $μBeff=μB−Δμ$⁠. In particular, Fig. 2(b) shows that increasing $μBeff$ in the CG model leads to a non-monotonic behavior in the variance of the particle number N (analogous to the compressibility of the model). This maximum indicates that the system has a tendency for de-mixing at larger $ηSr$ (one expects a divergent compressibility at the critical point, if one exists). In the following, we fix μ_B at the value corresponding to this maximum—the relatively large fluctuations at this point are challenging for the multi-level model because the distributions $pC(C)$ and $pF(C)$ are broader, requiring good sampling. The corresponding CG system has an average of N ≈ 11.6 big particles, occupying around 20% of the available volume.

For the specific quantities that we will compute for this mixture, Fig. 3 shows the expectations of the big-particle pair correlation function g(r) and the distribution of the number of big particles P(N). Results are shown for both CG and FG models (in the FG case, the results are computed using the two-level method). For both quantities of interest, the CG model provides an accurate but not exact description of the model. In particular, the CG model underestimates the pair correlation at the point where two big particles are in contact. The distributions of the number of big particles in Fig. 3(b) are both unimodal: both the FG and CG systems are well below the critical point of demixing.

Compared to the critical hard-sphere mixture discussed in Ref. 17, the system we consider here is smaller and has a lower volume fraction $ηSr$ of the small particles. This is still challenging for conventional Monte Carlo algorithms but can be simulated fast enough to evaluate the performance and compare the computational methods discussed here. Furthermore, the lower small-particle volume fraction helps with the construction of the intermediate level in Sec. IV, whose underlying approximation decays as $ηSr$ increases (see  Appendix B).

This section reviews the two-level method of Refs. 17 and 18 and then lays out its three-level extension. The presentation of the method is intended to be generic and applicable to a variety of systems. However, we first introduce the key ideas using the example and illustrations of Fig. 1 for the hard-sphere mixture.

The two-level method is constructed with the scale separation of the mixture in mind: it splits the simulation of the big and small spheres into two stages by first simulating a CG system of large particles alone and computing ⟨A⟩_C. Then, differences between ⟨A⟩_C and ⟨A⟩_F are computed by a reweighting (importance sampling) method. The weight factors for this computation are obtained by an annealing step, where the small particles are slowly inserted into the system, with the large particles held fixed [see Fig. 1(e)]. The advantage of this procedure is that large particle motion only happens in the CG simulation where the small particles are absent—there is no scale separation in this case so simulations are tractable. Similarly, insertion of the small particles happens in a background of fixed large particles, so these annealing simulations do not suffer long time scales associated with large-particle motion. This makes for tractable simulations in scale-separated systems, as long as the CG model is sufficiently accurate: see Refs. 17 and 18 for further discussion.

In practice, the simulation effort for two-level computations is dominated by the annealing step. The weighting factors are required to high accuracy, which means that the annealing must be done gradually. Moreover, the weights are subject to numerical uncertainties that tend to be large in systems with many small particles. This limits the method to systems of moderate size, with moderate $ηSr$ (see Ref. 17).

We show in this work that such problems can be reduced by breaking the annealing process into several stages—this is the idea of the three-level method [Fig. 1(e)]. Specifically, we start (as before) with a population of configurations of the CG model. We perform a first annealing step where the small particles are added in regions that are close to large ones. The information from this step is used in a resampling process, which partially corrects the coarse-graining error by discarding some of the configurations from the population and duplicating others. (This idea is similar to go-with-the-winners.²⁷) Finally, the second annealing step inserts the small particles in the remaining empty regions, arriving at configurations of the FG model. Hence, the end point is the same as the two-level method, but the annealing route is different.

In practice, the effectiveness of the three-level method relies on a clear physical understanding of the intermediate (partially inserted) system in order to decide which configurations to discard in the resampling step. For the hard-sphere case, that issue will be discussed in Sec. IV; a more general discussion is given in Sec. VII. The remainder of this section describes the two- and three-level methods in more detail.

We review the two-level method of Refs. 17 and 18. For a general presentation, we assume that CG and FG models exist with configurations $C$ and $X=(C,F)$⁠, respectively. In the case of hard spheres, $C$ and $F$ correspond to configurations of the large and small spheres, respectively.

The two-level method is an importance sampling⁴¹ (or reweighting) computation, closely related to the free-energy perturbation method of Zwanzig.¹⁹ We use the grand canonical Monte Carlo (GCMC) method to sample M_C configurations from p_C, and these are denoted by $C1,C2,…,CMC$⁠. Then, the CG average can be estimated as

As the sampling is increased (M_C → ∞), we have $ÂC→⟨A⟩C$⁠. However, if the coarse-graining error

is significant, $ÂC$ does not provide an accurate estimate of ⟨A⟩_F.

To address this problem, we use an annealing procedure based on Jarzynski’s equality²⁰ that starts from a coarse configuration $C$ and populates the fine degrees of freedom $F̂$⁠; at the same time, it generates a random weight $Ŵ(C)$ with the property that

where the angle brackets with subscript J indicate an averaging over the annealing process (analogous to Jarzynski’s equality²⁰), and ξ is a constant (independent of $C$⁠). The details of the annealing process are given in  Appendix A. It is applied to a set of M_F coarse configurations, again denoted by $C1,C2,…,CMF$⁠, which are typically a subset of the M_C CG configurations above.

For later convenience, we define

In practical applications, the constant ξ is not known, but its effect can be controlled by defining the self-normalized weight

Since the $Cj$ are representative of p_C, the denominator in $ŵ$ converges to ξ as M_F → ∞ and so $ŵ(Cj)→Ŵn(Cj)$⁠. Then, the estimator

converges to ⟨A⟩_F as M_F → ∞. (In the case that $Ŵ$ is not random then this procedure recovers the free energy perturbation theory of Zwanzig.¹⁹)

The annealing process has one useful additional property: Let the joint probability density for the weight and the fine degrees of freedom be $κ(Ŵ,F̂|C)$⁠, which is normalized as $∫κ(Ŵ,F̂|C)dF̂dŴ=1$⁠. We show in  Appendix A that

This formula is the essential property of the annealing procedure, which is required for the operation of the method. Additionally integrating over $F̂$ shows that (14) ensures that (10) also holds. This means in turn that if $B=B(C,F)$ is an observable quantity that depends on both coarse and fine degrees of freedom,

converges to ⟨B⟩_F as M_F → ∞.

This method can be easily improved without extra computational effort. The key idea^31–33 is to estimate the FG average as the sum of the CG average and the coarse-graining error (9),

Then, use importance sampling to estimate Δ, as

Finally, a suitable estimator for the FG average is obtained by combining the estimate of the coarse-graining error with the corresponding CG quantity,

This estimator converges to ⟨A⟩_F in the limit where M_C, M_F → ∞. As discussed in Ref. 18, the variance of the estimate $Δ̂$ is typically smaller than that of $ÂF$⁠, and the CG estimate $ÂC$ is cheap to compute accurately. Thus, the combined difference estimator $ÂF,Δ$ is typically more accurate at fixed computational cost.

The importance sampling methodology has a useful physical interpretation, which we explain for the example of the hard-sphere mixture. If we consider a fixed configuration of the large particles, the grand canonical partition function for the small particles is

As the system is annealed (the small particles are inserted), we estimate (19) by a free-energy method based on Jarzynski’s equality²⁰ (see  Appendix A for details). Since the annealing is stochastic, this yields an estimate of the partition function, which we denote by $Ξ̂[C,μS]$⁠. Moreover, this estimate is unbiased $⟨Ξ̂[C,μS]⟩J=Ξ[C,μS]$⁠. Hence, we can take

and using (4) and (5), we see that (10) holds, with ξ = (Ξ_C/Ξ_F).

Physically, the CG model is constructed so that the Boltzmann factor $e−UC(C)$ is a good estimate of the small-particle partition function $Ξ[C,μS]$⁠. If this is the case, the model is accurate. The two-level methodology uses estimates of the small-particle partition function (or, equivalently, their free energy) and compares it with the assumptions that were made about this quantity in the CG model. By analyzing the differences between these quantities, the differences between CG and FG models can be quantified. The effectiveness of this method for numerical simulation of the mixtures of large and small particles was discussed in Refs. 17 and 18.

The distribution of the importance weights $ŵ(C)$ impacts the accuracy of the resulting FG estimate $ÂF$⁠. Additionally, it serves as a useful indicator of the accuracy of the CG model and the variance of the free energy computation. To give an example, we apply the two-level method to the example problem from Sec. II C. In Fig. 4, we show the empirical distribution of 18 000 weights of the example system that are computed using an accurate annealing process; we use these computations as the reference solution in Sec. V. This illustrates a situation where the two-level method is applicable, where no single sample dominates and only very few samples have a weight larger than 10.

If one considers less accurate CG models, the variance of the weights increases, and the tail of their distribution gets heavier. Eventually, one would reach a situation where a few samples dominate the weighted sum (13). For accurately computed weights $ŵ(C)$⁠, such a breakdown of the two-level method indicates that the CG model is not sufficiently accurate. This behavior provides a useful feedback loop that can be used to iterate on the CG model itself.¹⁸ 

We now present the three-level method for estimation of $⟨A(C)⟩F$⁠.

We start by generating M₀ samples of the CG model, denoted by $C01,…,C0M0$⁠. The subscript indicates the step within the algorithm (which is 0 for the initial sampling of coarse configurations). The CG average of A can be estimated similarly to (8)

In addition to the CG and FG models, the three-level method also relies on an intermediate set of configurations, which correspond in the hard-sphere mixture to the system where the small particles have been inserted in regions close to the large ones (see Fig. 5). This state is described by an equilibrium probability distribution

where $UI(C,F)$ is an interaction energy. Its construction for the hard-sphere mixture will be discussed in Sec. IV.

The first annealing step of the three-level algorithm applies the two-level method, with the FG distribution p_F replaced by p_I. This part of the algorithm closely follows Sec. III B, and we give a brief discussion that mostly serves to fix notation. We start with a set of M₁ coarse configurations that are samples of p_C; they are denoted by $C11,C12,…,C1M1$⁠, where the subscript 1 indicates the intermediate stage of the three-level method. (These will typically be a subset of the configurations that were generated on the coarse level.)

For each coarse configuration $C1$⁠, we anneal the fine degrees of freedom of the system $F̂1$ to arrive at the intermediate level and generate a random weight $Ŵ1(C1)$ with the property

with a constant ξ₁ independent of $C1$⁠. (For the hard spheres, we recall that particles are inserted preferentially in regions close to large ones, and this is illustrated in the top row of Fig. 5.)

As before, we define $Ŵ1n(C1)=Ŵ1(C1)/ξ1.$ Again the constant ξ₁ is generally not known, so we define the self-normalized weight

which converges to $Ŵ1n(C1j)$ as M₁ → ∞. Then, the estimator

converges to ⟨A⟩_I as M₁ → ∞. Similar to (14), the joint probability density $κ1(Ŵ1,F̂1∣C1)$ of the weight and fine degrees of freedom at the intermediate level, defined by the annealing process, fulfills

Hence, similar to (15), we also obtain

which converges to ⟨B⟩_I as M₁ → ∞.

At the end of the intermediate level, we have M₁ large-particle configurations. For each configuration $C1j$⁠, the process of annealing to the intermediate level also provided the weight $ŵ1(C1j)$ and the small-particle configuration $F̂1j$⁠. This information can be used to build a set of configurations that are representative of p_I. This procedure is called resampling, its validity in this example relies on the property (26) of the annealing procedure. This is the part of the method that is similar to population-based sampling approaches, such as SMC²⁸ or go-with-the-winners.²⁷ The idea is that one should focus the effort of the annealing process onto coarse configurations that are typical of the full system, and to discard those which are atypical, see Fig. 6 for a visualization of this step.

We write $X̂1j=(C1j,F̂1j)$ for the full configuration that is obtained by the annealing procedure at the intermediate level. The resampled configurations will be denoted by $X21,X22,…,X2M2$⁠; they are representative of the intermediate level p_I. There are M₂ of them, and the subscript 2 indicates the final stage of the three-level method. The simplest resampling method (multinomial resampling) is that each $X2i$ is obtained by copying one of the $X̂1j$⁠, chosen at random with probability $ŵ1(C1j)$⁠. In applications, one typically replaces this by a lower variance resampling scheme like residual resampling; see Ref. 42 for a comparison of commonly used variants.

We then perform the second annealing step that starts from an intermediate level configuration $X2=(C2,F2)$ and anneals the fine degrees of freedom from the intermediate to the fine level, yielding $F̂2$ and a weight $Ŵ2(X2)$⁠; details are given in  Appendix A. For the hard sphere system, this involves further insertion of small particles to fill the system and generate realistic configurations of the full mixture. This procedure is shown in the bottom row of Fig. 5.

Since the starting point of the annealing procedure is $X2$⁠, the joint probability density of the annealing process $κ2(Ŵ2,F̂2∣X2)$ depends on both large and small particles. Therefore, the analog of (26) requires an additional average over the small particles of the starting configuration

for some constant ξ₂. Note that $pI(F|C)=pI(C,F)/pI(C)$⁠. Similar to (23), the weights $Ŵ2(X2)$ have the property

From here, we proceed as before. We define the normalized weight $Ŵ2n(X2)=Ŵ2(X2)/ξ2$ and its self-normalized estimate

Since the $X2j$ are representative of p_I, it follows from (28) that observables of the coarse system A can be estimated as

which converges to ⟨A⟩_F as M₂ → ∞. Similar to (15), we can also obtain a consistent FG estimates of observable quantities B that depend both on coarse and fine degrees of freedom by

Following the same variance reduction strategy as in Sec. III B, we can define a difference estimator of the FG average, which is expected to have lower statistical uncertainty: let

Then,

is a consistent estimator of ⟨A⟩_F, analogous to (18).

A few comments on the three-level method are in order. First, there is a simple generalization to four or more levels by splitting the annealing procedure into more than two stages. As such, the method is an example of a sequential Monte Carlo (SMC) algorithm (which is sometimes more descriptively referred to as sequential importance sampling and resampling^28,43–45). We note from (10) that the weights obtained from the annealing step are random, and this is not the standard situation in SMC, but similar ideas have been previously studied in Refs. 46–49. Combining an SMC algorithm with a difference estimate as in (35) has been investigated in Refs. 34–36.

Second, we observe that the key distinction between the two- and three-level algorithms is the resampling step at the intermediate level. Without this, the three-level method reduces to a simple two-level method with an arbitrary stop in the middle of the annealing process. As noted above, the resampling process is designed to partially correct differences between the CG and FG models. This relies on a good accuracy of the intermediate level (otherwise the wrong configurations might be discarded, which hinders numerical accuracy). On the other hand, we note that for sufficiently large numbers of samples M₀, M₁, M₂, the method does provide accurate FG estimates, even if the CG and intermediate level models are not extremely accurate. The distinction between the different methods comes through the number of samples that are required to obtain accurate FG results.

Third, note that the ideal situation for difference estimation is that the three terms in (35) get successively smaller. That is, the coarse estimate is already close to ⟨A⟩_F, the intermediate-level estimate provides a large part of the correction, and the fine-level correction is small. In this case, it is natural to use a tapering strategy where the number of samples used at each level decreases,

This allows a fixed computational budget to be distributed evenly between the various levels to minimize the total error.

As noted above, the intermediate probability distribution p_I must be designed carefully in order for the resampling part of the three-level method to be effective. We now describe how this is achieved for the hard sphere mixture.

To motivate the intermediate level, recall Fig. 5, and note that defining a suitable CG model is equivalent to an estimate of the small-particle free energy in the final (fully inserted) state. The physical idea of the intermediate level is that the free energy associated with the first stage of insertion may be hard to estimate (because of the complicated packing of the small particles around the large ones), but the free energy difference associated with the second stage should be easier (because it corresponds to insertion into large empty regions where the packing of the small particles is similar to that of a homogeneous fluid, whose free energy can be estimated based on analytic approximations). A combination of these ideas yields an intermediate level that represents the big particle statistics more accurately than the CG model. Similar ideas have been considered before in multi-scale simulation,^4,50,51 in particular the problem of estimating the small-particle free energy has some similarities to estimation of solvation free energies (where the depletant here plays the role of a solvent).

We start by analyzing the small particles, so we fix the large particles in some configuration $C$⁠. The idea of the intermediate level is to first insert small particles only in a region close to the large particles $C$⁠, and then use this information to make the intermediate marginal distribution $pI(C)$ match the FG marginal $pF(C)$ as closely as possible. The structure of the intermediate level is depicted in the bottom row of Fig. 1(e), and an example configuration is shown in Fig. 1(d). We implement this idea by introducing an effective (one-body) potential that acts on the small particles. We first define

to be the distance from the point r to the nearest large particle. Small-particle insertion is suppressed in regions far from large particles by a potential energy term

where

and the function

interpolates from zero (for small distances r) to the value s at large r. This function acts as a smoothed out step function, where δ_free is the position of the step and l its width. In Figs. 1(e) and 5, areas where $EC(r)>0$ are indicated by blue shaded regions, in which the insertion of small particles is suppressed.

Then, define a grand-canonical probability distribution for the small particles in the partially inserted (intermediate) system as

This distribution is normalized as $∫p̃I(F∣C)dF=1$⁠. It depends on the three parameters s, δ_free, l, as well as the underlying parameters of the hard sphere mixture model.

The next step is to construct the weights $Ŵ1(C)$⁠. For consistency with (22), we write the intermediate-level distribution in the form

As discussed above, the term Φ^corr should be designed so that the respective coarse-particle marginals $pI(C)$ and $pF(C)$ match as closely as possible. Using (4), (19), and (41), we can show that a perfect match requires $Φcorr(C)=Φex(C)$ with

where ϕ₀ is an irrelevant constant. Since the Ξs in (43) are partition functions, determination of Φ^ex reduces to computation of the free energy difference between the non-homogeneous small particle distributions of the partially and fully inserted system. We now explain how Φ^corr is defined, as an approximation to Φ^ex.

As a preliminary step for estimating Φ^ex, we first consider the grand potential Φ for the small particles, in a system with no large particles, where the small particles feel an (arbitrary) smooth potential $E=E(r)$⁠. The grand potential of this system is

where 0 indicates the large-particle configuration with no particles at all (N = 0).

If $E$ varies slowly in space, a simple approach to this integral is to assume that the system is locally the same as a homogeneous system in equilibrium—similar to the local density approximation.⁵² In this case,

where $p$ is the pressure, expressed as a function of the chemical potential.

However, this approximation is not sufficiently accurate for the current application. To this end, we include a correction to account for inhomogeneities, as a squared gradient term $Φ[E;μS]≈Φsq[E;μS]$ with

(Within a gradient expansion, this is the first correction that is consistent with rotational and inversion symmetry.)

We show in Appendix B 1 that $g$ can be estimated as

where S(μ; q) is the structure factor of the small hard-sphere system. For a numerical estimate of this Φ^sq, we estimate the pressure $p$ by the accurate equation of state from Ref. 38, and $g$ is estimated from (47) using the structure factor from Ref. 53. A numerical example demonstrating the accuracy of this second order approximation for a non-homogeneous hard-sphere fluid can be found in Appendix B 2.

We are now in a position to approximate Φ^ex in terms of Φ^sq. This (analytical) calculation is illustrated in Fig. 7. We require an estimate of Φ^ex, which is the free-energy difference between the partially inserted and fully inserted systems in panels (a) and (b). This is achieved as a sum of three free-energy differences. In the first step, the large particles are removed and the small-particle fluid is re-equilibrated, to fill up the remaining space, leading to panel (c). Then, the confining potential $Ũ$ is removed and the small particles fully inserted, leading to (d). Finally, the large particles are re-inserted and the small particles re-equilibrated again, leading to (b).

To make this precise, define $ΦC[E;μS]$ as the grand potential of the small particles in the potential $E$⁠, where the large particles are also included, with configuration $C$⁠. Then, the desired free energy difference between panels (a) and (b) is

where we took ϕ₀ = 0.

From the definitions in Sec. IV A, the free energy difference between panels (c) and (d) is $Φ[0;μS]−Φ[EC;μS]$⁠, from (38) and (44). Our central approximation is that the free energy difference between panels (a) and (c) is (approximately) equal and opposite of the difference between (d) and (b) because the local environment of the large particles is the same in both cases. (The only differences are in regions far from any large particles.) At this level of approximation, the free energy differences between [(a),(b)] and [(c),(d)] are equal,

Finally, the right-hand side can be estimated by the square gradient approximation (46), yielding $Φex(C)≈Φcorr(C)$ with

Operation of the three-level method requires numerical estimates of this Φ^corr, which includes the integral in (46). Moreover, its value is exponentiated when computing weight factors $Ŵ$⁠, so these numerical estimates are required to high accuracy. This is a non-trivial requirement because the integrand is constant on regions far from the big particles, but it varies much more rapidly when these particles are approached. In such situations, adaptive quadrature schemes are appropriate: we use the cuhre algorithm of the cuba library⁵⁴ that uses globally adaptive subdivision to refine its approximations in the relevant regions of space. Note, however, that while the choice of the numerical integrator influences the intermediate level, small errors in estimation of this integral will be corrected by the second annealing step, so such errors do not affect the consistency of our numerical estimators.

Given this choice of $Φcorr(C)$⁠, the intermediate level distribution p_I of (42) has been completely defined, although it still depends on the three parameters δ_free, s, l that appear in the function ɛ(r). We also note that given the approximations made, it is not expected that this p_I is optimal (its marginal $pI(C)$ does not match $pF(C)$ perfectly). Subsection IV C discusses the parameter choices, and some possibilities for correction factors that can be added to Φ^corr in order to address specific sources of error.

In fixing the parameters δ_free, s, l, several considerations are relevant. First, if s is too small or δ_free is too large, the potential $EC$ has little effect on the system and the small particles are not restricted to be close to the large ones. In this case, p_I ends up close to p_F, and there is little benefit from the intermediate level. On the other hand, the accuracy of Φ^sq is greatest when the gradient of the potential $EC$ is small, and this favors small s and large l, δ_free. In practice, it is also convenient if the two annealing stages insert similar numbers of particles, so that their computational costs are similar. For the example system in Sec. II C, we will present results for a suitable parameter set

We have also tested other values, a few comments are given below.

We will consider several variants of the intermediate level. We denote by $pI(1)$ the distribution defined by (42) and (50), with parameters (51). Figure 8 shows how the quantities of interest differ between the CG and FG models, and the corresponding differences between the intermediate level and the FG model. Here, Δg(r) is the difference between g(r) for the FG model and the distribution of interest (which is either the CG distribution p_C or one of the variants of the intermediate distribution). ΔP(N) is the corresponding difference in the probability that the system has N large particles.

For the value of g(r) at contact, we see that the intermediate level $pI(1)$ corrects around half of the deviation between CG and FG models. However, the probability distribution of the N has the opposite situation that the intermediate level is less accurate than the CG model. [This is partly attributable to the fact that Δμ in Eq. (6) has been chosen to make the CG model accurate.]

To explore the behavior of the intermediate level, we constructed two variants of p_I. The aim is to understand why $pI(1)$ has inaccuracies and to (partially) correct for them. There are two main approximations in the intermediate level $pI(1)$⁠: the first is (49) and the second is that $Φ[EC,μS]$ can be approximated by the square-gradient approximation (46). The first approximation neglects a significant physical phenomenon in these systems, which is a layering effect of the small particles around the large ones. This is illustrated in Fig. 9 by the radial distribution function $gBS0$ between large and small particles (measured in a system with a single large particle). One sees that there is typically an excess of small particles close to the large ones, followed by a deficit $[gBS0(r)<1]$⁠, and a (weak) second layer.

For (49) to be accurate, the intermediate level should have enough small particles to capture this layering so that the particles being inserted in the second annealing stage are not strongly affected by the presence of the large particles. However, computational efficiency requires that δ_free is not too large, so these layers are not fully resolved at the intermediate level. To partially account for this effect, we make an ad hoc replacement of μ_S in (46) by an effective chemical potential μ_lay(r), which is chosen such that the corresponding reservoir volume fraction $ηSlay(r)$ satisfies

In estimating the free energy of the small particles that are inserted in the second level of annealing, this adjustment to Φ^sq helps to counteract the error made in (49), leading to an updated potential Φ^corr,2. The intermediate level constructed in this way is denoted by $pI(2)$⁠. The results of Fig. 8 show that this variant is (somewhat) more accurate than $pI(1)$⁠. However, the intermediate level still tends to have a smaller number of large particles than the full (FG) mixture.

To investigate this further, we took 800 representative CG configurations. For each one, we estimate the error associated with the approximation (50)

Results are shown in Fig. 10. One sees that the errors are of order unity (note that Φ^ex itself is of order 10⁴ so this is a small relative error; see below); there is a systematic trend that Φ^corr,2 underestimates Φ^ex when N is large. To correct this error, we introduce an additional correction term to Φ^corr,2,

and denote the intermediate level constructed in this way by $pI(3)$⁠.

A least squares fit to Fig. 10 suggests to take α_corr = 0.076; in practice, this tends to over-correct the error in $Φcorr,2(C)$⁠, and we find better performance with a smaller value

However, the performance of the method depends only weakly on the specific choice of α_corr, which is discussed in  Appendix C. For all following results, we define the intermediate level $pI=pI(3)$ to use the potential Φ^corr,3.

An important aspect of the three-level method is the self-consistency of the general approach. The intermediate level variants $pI(1)$ and $pI(2)$ were constructed on a purely theoretical basis. The corresponding results in Fig. 8 indicated good performance, but that the distribution of N had a systematic error. This error was quantified precisely in Fig. 10, which enabled an improvement to the intermediate level. In principle, this procedure could be repeated to develop increasingly accurate variants of p_I. That approach would be useful if (for example) one wanted to consider increasingly large systems, where the requirements for the accuracy of p_I become increasingly demanding.

One way to see the effect of system size is to note that Fig. 10 required the estimation of $ΦC[EC,μS]$ and $ΦC[0,μS]$⁠, whose values are of order 1 × 10⁴. Since the free energies are exponentiated in the weights for resampling, an absolute error of ±1 is required on these free energies, while their absolute values are extensive in the system size. Hence, one sees that accurate estimates of the free energy are required: their relative error is required to be of the order of the inverse volume of the system.

In this section, we apply the three-level method to the example from Sec. II C using the intermediate level from Sec. IV, with the parameters defined in (51) and (55). The parameters and the annealing schedules are chosen such that, on average, the first and second steps have the same computational effort; see  Appendix A for details.

It can be proven⁴⁹ that the three-level method provides accurate results, in the limit where the population sizes M₀, M₁, M₂ are all large. In particular, we expect the estimators $ÂF3L,ÂF,Δ3L$ to all obey central limit theorems (CLTs), and the two-level estimators $ÂF,ÂF,Δ$ behave similarly. Detailed results are given in Sec. VI. The important fact is that, for large populations, the variances of the estimators behave as

where M is the relevant population size and Σ is called the asymptotic variance (it depends on the observable A and on which specific estimator is used). In general, the estimators may have a bias, which is also of order 1/M. This means that the uncertainty in our numerical computations is dominated by the random error, whose typical size is $Σ/M$⁠, and the mean squared error is given by the variance $MSE(Â)=Var(Â)$⁠, to leading order.

This gives us an easy way to measure and compare the performance of the different estimators. Suppose that we require an estimate of A with a prescribed mean squared error. The associated computational cost can be identified with the population size M and is given by (56) as $M≈Σ/MSE(Â)$⁠. Clearly, estimators with small Σ should be preferred. In practice, we do not compare computational costs at fixed error, instead we compare variances $Var(Â)$ at fixed M. For any two algorithms (and assuming that M is large), the ratio of these variances approximates the ratio of the Σ’s, which can then be interpreted as a ratio of computational costs (at fixed MSE). Numerical results are presented in Sec. V B.

We note that the theoretical results for convergence do not require that the coarse or intermediate levels are accurate. However, one easily sees¹⁸ that serious inaccuracies in these levels lead to very large Σ. In such cases, one may require prohibitively large populations to obtain accurate results.

In this section, we demonstrate (for the example in Sec. II C) that we do not require very large populations for the three-level method and that the numerical results are consistent with (56). After that, we estimate the asymptotic variances for the two-level and three-level methods. We will find that introducing the third level improves the numerical performance, corresponding to a reduction in Σ.

To this end, we investigate the pair correlation g(r) of the big particles. As seen in Fig. 8(a), the coarse approximation of g(r) has a substantial error, especially when two big particles are in contact. To quantify this specific effect, we define the coordination number N_c, which is the number of large particles within a distance r₁ of a given large particle. (For a given configuration, this quantity is estimated as an average over the large particles. We take r₁ ≈ 10.73σ_S to be the first minimum of g(r) of the CG model.) For our example, the coordination number for the FG and CG systems are given by

To illustrate the reliable performance of the method, we take a simple example with M₀ = 4 × 10⁵ and M₁ = M₂ (no tapering) and we focus on the difference estimator $ÂF,Δ3L$⁠, which we expect to be the most accurate. The corresponding numerical estimate of g(r) is denoted by $ĝ(r)$⁠, binned using 40 equidistant bins at positions r_j between r = 10 and r = 12. Figure 11(a) shows estimates of the difference between $ĝ(r)$ and its true value, as the population size increases. (The FG result was estimated independently by the two-level method using a large value of M_F = 18 000.) A population M₂ of several thousand is sufficient for an accuracy better than 0.5 in each bin of g(r).

For smaller M₂, fluctuations in the measured $ĝ(r)$ are apparent in Fig. 11(a). To estimate their size, we define the error for a single run of the three-level method by summing over the bins,

Hence, one expects from (56) that this error decays with increasing population, proportional to $MF−1/2$⁠. Figure 11(b) shows an estimate of (58), which is consistent with this expected scaling.

We now investigate whether the three-level method does indeed improve on the performance of the (simpler) two-level method of Refs. 17 and 18. The key question is whether the resampling step is effective in focusing the computational effort on the most important configurations of the big particles.

We recall from above that removing the resampling step from the three-level method leads to a two-level method, where the annealing process is paused at the intermediate level and then restarted again. In order to test the effect of resampling, we compare these two schemes, keeping the other properties of the algorithm constant, including the annealing schedule. (To test the overall performance, one might also optimize separately the annealing schedules for the two-level and three-level algorithms and compare the total computational time for the two methods to obtain a result of fixed accuracy. However, such an optimization would be very challenging, so instead we focus on the role of resampling.)

As a very simple quantity of interest, we take the co-ordination number N_c. We run the whole algorithm N_runs independent times, and we estimate N_c for each run. This can be done using several different estimates of N_c. These are (i) the two-level estimates $ÂF$ and $ÂF,Δ$ from (13) and (18); (ii) the corresponding three-level estimates $ÂF3L$ and $ÂF,Δ3L$ of (31) and (35), in which we also vary the ratio M₁: M₂, to see the effects of tapering.

All comparisons are done with a fixed total computational budget. We have chosen parameters such that the first and second annealing stages have the same (average) computational cost. This means we need to hold M_T = (M₁ + M₂)/2 constant during tapering. The two-level method takes M_F = M_T (because the single step of annealing in the two-level method has the same cost as the two annealing steps of the three-level method). For the coarse level estimates $ÂC$ and $ÂC3L$ (which are used in computation of $ÂF,Δ$ and $ÂF,Δ3L$⁠), the CG computations are cheap so we take M₀ = M_C = 6 × 10⁶. This is large enough that the numerical errors on these coarse estimates are negligible in comparison to the errors from higher levels.

For each version of the algorithm, we measure the sample variance of the N_runs estimates. Results are shown in Fig. 12 for N_runs = 60 and M_T = 500. The error bars are computed by the bootstrap method.⁵⁵ It is useful that the variance of all these estimators is expected to be proportional 1/M_T: this means that reducing the variance by a factor of α requires that the computational effort is increased by the same factor. Hence, the ratio of variances of two estimators is a suitable estimate for the ratio of their computational costs.

When carrying out these runs, each estimator was computed by performing annealing on the same set of coarse configurations, to ensure direct comparability. (More precisely, we take a set of 700 representative configurations that are used for the method with M₁:M₂ = 7:3, and other versions of the method used a subset of these 700.) In addition, it is possible to share some of the annealing runs when computing the different estimators (while always keeping the 60 different runs completely independent). This freedom was exploited as far as possible, which reduces the total computational effort. However, it does mean that the calculations of the different estimators are not at all independent of each other.

All three-level estimators have a reduced standard deviation compared to their two-level equivalents, demonstrating the usefulness of the intermediate resampling step. In all cases, the difference estimate outperforms its equivalent final-level estimate; this effect is stronger for the three-level estimate, providing evidence that the intermediate stop additionally improves the quality of the control variate in the difference estimates.

The effect of introducing tapering from M_I = 600 to M_F = 400 is difficult to assess, given the statistical uncertainties in this example. The variance of the tapered final-level estimator is very close to the non-tapered one, despite averaging over fewer configurations. This is possible since we start with more samples in the CG model that improves the sampling at the intermediate step, where we then resample to keep relevant configurations. As the results for the 700–300 tapering show, the tapering rate needs to be chosen carefully as a too aggressive rate can decrease the performance quickly.

Overall, the numerical tests in this section provide strong evidence of the benefit of the intermediate resampling. For our example, switching from a two-level to a three-level difference estimator substantially reduces the variance, from around 0.0029 for the two-level method to 0.0016 at a fixed computational budget. As discussed just below (56), the ratio of these numbers can be interpreted as the ratio of costs for the two- and three-level method: the conclusion for this case is that including the intermediate level reduces the cost by ∼45%. This demonstrates a significant speedup in this specific case, which provides a proof-of-principle of the approach.

In Sec. V, we have seen that the three-level method outperforms the two-level method in numerical tests, both for the final-level as well as the difference version of the estimator. In this section, we provide convergence results for both algorithms and compare their asymptotic performance as the number of configurations goes to infinity.

The proof is general, but it does require some assumptions on the models of interest. First, for every allowed CG configuration [that is, configurations $C$ with $pC(C)>0$], we assume that the quantity of interest A is bounded. In addition, the probability density $pC(C)$ must be non-zero whenever $pI(C)$ is non-zero, and similarly $pI(C,F)$ must be non-zero whenever $pF(C,F)$ is non-zero.

The two-level method has been previously analyzed in Ref. 18. We summarize its key properties. It was noted in Sec. III B that $ÂF→⟨A⟩F$ as M_F → ∞ (specifically, this is convergence in probability⁵⁶). We also expect a CLT for this quantity: as in Eq. (56), the distribution of the error $(ÂF−⟨A⟩F)$ converges to a Gaussian with mean zero and variance Σ_F/M_F. We will derive a formula for this variance, which will be compared later with the corresponding quantity for the three-level model.

For compact notation, it is convenient to define the re-centered quantity of interest