Correction of coarse-graining errors by a two-level method: application to the Asakura-Oosawa model

TL;DR

This paper introduces a two-level correction method that combines coarse-grained and fine-grained simulations to improve the accuracy of physical property predictions, demonstrated on the AO colloid-polymer model.

Contribution

The paper presents a novel self-consistent two-level simulation approach to correct coarse-graining errors in physical system models.

Findings

The liquid-vapour critical point is influenced by three-body interactions.

The method's accuracy depends on the computational effort.

Three-body interactions are significant in the AO model.

Abstract

We present a method that exploits self-consistent simulation of coarse-grained and fine-grained models, in order to analyse properties of physical systems. The method uses the coarse-grained model to obtain a first estimate of the quantity of interest, before computing a correction by analysing properties of the fine system. We illustrate the method by applying it to the Asakura-Oosawa (AO) model of colloid-polymer mixtures. We show that the liquid-vapour critical point in that system is affected by three-body interactions which are neglected in the corresponding coarse-grained model. We analyse the size of this effect and the nature of the three-body interactions. We also analyse the accuracy of the method, as a function of the associated computational effort.

Equations

$$p_{\mathrm{f}}(\mathcal{C},\mathcal{F})=\frac{1}{Z}{\mathrm{e}}^{-\beta E(\mathcal{C},\mathcal{F})}$$

$$⟨A{⟩}_{\mathrm{f}}=\int p_{\mathrm{f}}(\mathcal{C},\mathcal{F})A(\mathcal{C})\mathrm{d}\mathcal{C}\mathrm{d}\mathcal{F}.$$

$$p_{\mathrm{m}}(\mathcal{C})=\int p_{\mathrm{f}}(\mathcal{C},\mathcal{F})\mathrm{d}\mathcal{F}.$$

$$p_{\mathrm{c}}(\mathcal{C})=\frac{1}{Z_{\mathrm{c}}}{\mathrm{e}}^{-\beta E_{\mathrm{c}}(\mathcal{C})}$$

$$⟨A{⟩}_{\mathrm{c}}=\int p_{\mathrm{c}}(\mathcal{C})A(\mathcal{C})\mathrm{d}\mathcal{C}.$$

$${\hat{A}}_{\mathrm{c}}=\frac{1}{{\mathcal{N}}_{\mathrm{c}}}\sum _{i=1}^{{\mathcal{N}}_{\mathrm{c}}}A({\mathcal{C}}_i).$$

$$⟨A{⟩}_{\mathrm{f}}-⟨A{⟩}_{\mathrm{c}}$$

$$={⟨[w(\mathcal{C})-1]A(\mathcal{C})⟩}_{\mathrm{c}}$$

$$\hat{\Delta }=\frac{1}{{\mathcal{N}}_{\mathrm{f}}}\sum _{i=1}^{{\mathcal{N}}_{\mathrm{f}}}[\hat{w}({\mathcal{C}}_i)-1]A({\mathcal{C}}_i)$$

$$\hat{A}={\hat{A}}_{\mathrm{c}}+\hat{\Delta }.$$

$$\hat{A}\to ⟨A{⟩}_{\mathrm{f}}.$$

$$p_{\mathrm{f}}(\mathcal{C},\mathcal{F})=\frac{1}{{\Xi }_{\mathrm{f}}}{\mathrm{e}}^{\beta {\mu }_{\mathrm{B}}N+\beta {\mu }_{\mathrm{S}}n-\beta U_{\mathrm{f}}(\mathcal{C},\mathcal{F})}$$

$$q=({\sigma }_{\mathrm{S}}/{\sigma }_{\mathrm{B}}).$$

$$p_{\mathrm{m}}(\mathcal{C})=\frac{1}{{\Xi }_{\mathrm{f}}}\exp \left[\beta {\mu }_{\mathrm{B}}N+\frac{6{\eta }_{\mathrm{S}}^{\mathrm{r}}}{\pi {\sigma }_{\mathrm{S}}^3}{\mathcal{V}}_{\mathrm{a}}(\mathcal{C})\right],$$

$$\beta U_{\mathrm{c}}(\mathcal{C})=\sum _{i=1}^{N-1}\sum _{j=i+1}^N\beta V_{\mathrm{AO}}(|{\mathbf{R}}_i-{\mathbf{R}}_j|)$$

$$\beta V_{\mathrm{AO}}(r)=\frac{-{\eta }_{\mathrm{S}}^{\mathrm{r}}}{q^3}\left[(1+q{)}^3-\frac{3r(1+q{)}^2}{2{\sigma }_{\mathrm{B}}}+\frac{r^3}{2{\sigma }_{\mathrm{B}}^3}\right].$$

$$p_{\mathrm{c}}(\mathcal{C})=\frac{1}{{\Xi }_{\mathrm{c}}}{\mathrm{e}}^{\beta {\mu }_{\mathrm{B}}N-\beta U_{\mathrm{c}}(\mathcal{C})},$$

$$\frac{p_{\rm m}(\mathcal{C})}{p_{\rm c}(\mathcal{C})}\propto\exp\Big{(}\beta N\Delta\mu-\sum_{ijk}\beta V_{3}(\bm{R}_{i},\bm{R}_{j},\bm{R}_{k})\\ -\sum_{ijkl}\beta V_{4}(\bm{R}_{i},\bm{R}_{j},\bm{R}_{k},\bm{R}_{l})-\dots\Big{)}$$

$$P_{\mathrm{b}}(N_{\mathrm{b}})=⟨{\delta }_{N,N_{\mathrm{b}}}{⟩}_{\mathrm{f}}.$$

$$⟨\hat{W}(\mathcal{C}){⟩}_{\mathrm{J}}=\exp \left[\frac{6{\eta }_{\mathrm{S}}^{\mathrm{r}}}{\pi {\sigma }_{\mathrm{S}}^3}{\mathcal{V}}_{\mathrm{a}}(\mathcal{C})+\beta U_{\mathrm{c}}(\mathcal{C})\right].$$

$$\hat{w}({\mathcal{C}}_i)=\frac{\hat{W}({\mathcal{C}}_i)}{\frac{1}{{\mathcal{N}}_{\mathrm{f}}}{\sum }_{j=1}^{{\mathcal{N}}_{\mathrm{f}}}\hat{W}({\mathcal{C}}_j)}.$$

$$⟨\hat{w}(\mathcal{C}){⟩}_{\mathrm{J}}\approx \frac{p_{\mathrm{m}}(\mathcal{C})}{p_{\mathrm{c}}(\mathcal{C})},$$

$${\Xi }_{\mathrm{S}}({\mu }_{\mathrm{S}},\mathcal{C})=\exp \left[{\mathrm{e}}^{\beta {\mu }_{\mathrm{S}}}{\mathcal{V}}_{\mathrm{a}}(\mathcal{C})/{\sigma }_{\mathrm{S}}^3\right]=\exp \left[\frac{6{\eta }_{\mathrm{S}}^{\mathrm{r}}}{\pi {\sigma }_{\mathrm{S}}^3}{\mathcal{V}}_{\mathrm{a}}(\mathcal{C})\right].$$

$$\mathcal{I}(\mathcal{C})={\int }_0^{M}n(t)\frac{\mathrm{d}}{\mathrm{d}t}\mu (t)\mathrm{d}t.$$

$$\mathcal{I}(\mathcal{C})=\sum _{j=1}^K{n}_j\Delta {\mu }_j$$

$$\hat{W}(\mathcal{C})={\Xi }_{\mathrm{S}}({\mu }_0,\mathcal{C}){\mathrm{e}}^{\beta U_{\mathrm{c}}(\mathcal{C})}\frac{1}{m}\sum _{p=1}^m{\mathrm{e}}^{\beta {\mathcal{I}}_p(\mathcal{C})}$$

$$\delta N=N-⟨N⟩,X=\frac{\delta N}{{⟨\delta N^2⟩}^{1/2}}.$$

$$\beta U_{2,3}(\mathcal{C})=\beta U_{\mathrm{c}}(\mathcal{C})+\beta U_3(\mathcal{C})$$

$$U_3(\mathcal{C})=\sum _{i=1}^{N-2}\sum _{j=i+1}^{N-1}\sum _{k=j+1}^N{V}_3(R_{ij},R_{ik},R_{jk})$$

$$\big{\langle}\hat{A}-\langle A\rangle_{\rm f}\big{\rangle}\to 0$$

$$h^2=⟨{\left(\hat{A}-⟨A{⟩}_{\mathrm{f}}\right)}^2⟩.$$

Full text

I Introduction

Coarse-grained models are useful throughout many areas of physics and chemistry Likos (2001); Pak and Voth (2018). Depending on their context, they are used in different ways. In some cases, the idea is to develop a simple model that captures important qualitative features of some physical system Hagan and Chandler (2006); Haxton and Whitelam (2012); Saric et al. (2014). Alternatively, one may develop a coarse-grained (CG) model that captures quantitatively the behaviour of a more detailed (fine-grained, FG) model Noid et al. (2008); Praprotnik, Delle Site, and Kremer (2007); Ouldridge, Louis, and Doye (2011); Mladek et al. (2013); Pak and Voth (2018). This second type of application is one of the foundations of multi-scale modelling approaches Karplus (2014); Warshel (2014).

Enormous theoretical and computational effort has been invested in the development of accurate CG models. In the simplest case, one starts from a FG model with a classical Hamiltonian, and one aims to characterise its potential energy surface (PES) as a function of some simpler collective variables Noid et al. (2008); John and Csányi (2017). Alternatively, the FG model might be a quantum mechanical system, and one develops a CG model that aims to fit the Born-Oppenheimer PES for the nuclear co-ordinates Karplus (2014); Warshel (2014). A number of recent works have used machine-learning methods Behler and Parrinello (2007); Bartók et al. (2010); Cheng et al. (2019) to derive approximate PESs in this case.

Independent of the specific application, the most common approach in such cases is to derive a CG model that is as accurate as possible, and then to estimate observable properties of the original (FG) model by computing the corresponding quantities in the CG model. If the CG model is an accurate representation of the fine one then this approach works well. However, it is very rare for the CG model to be perfect, which means that this approach almost always results in some systematic errors, when estimating properties of the fine model. No information about these errors is available from analysis of the CG model alone.

In this work, we discuss a method for eliminating these coarse-graining errors. Specifically, we make quantitative estimates of observable properties of the FG model, which we decompose as the value of the corresponding quantity for the CG model, and an explicit correction term. We show how this correction can be estimated, in an efficient way. Our method is inspired by recent mathematical studies, particularly the multi-level approach Heinrich (1998); Giles (2008); Anderson and Higham (2012), which has been developed recently for applications in Bayesian statistics and uncertainty quantification Hoang, Schwab, and Stuart (2013); Dodwell et al. (2015); Beskos et al. (2017). (We note that the term “multi-level” has a specific mathematical meaning and should not be confused with the “multi-scale” terminology used in physics and chemistry Karplus (2014); Warshel (2014).) Our method corresponds to a multi-level implementation with only two levels, so we call it a two-level (TL) method. It shares many features with importance sampling methods, biased sampling techniques Frenkel and Smit (2001); Bruce and Wilding (2003), and free-energy perturbation methods Zwanzig (1954) (see for example Ref. Cheng et al., 2019 for a recent application of this method, in order to correct for coarse-graining errors). Our implementation of the TL method is informed by recent mathematical results and has been chosen to minimise the systematic errors (bias) in our estimates of properties of the FG system. As we discuss below, we minimise the bias by using a method for free-energy estimation Hummer and Szabo (2001) that is based on the fluctuation theorems of Jarzynski Jarzynski (1997) and Crooks Crooks (2000). This free-energy estimation method is equivalent to the mathematical technique of annealed importance sampling Neal (2001).

The method that we present is general. As a specific application, we use it to analyse the Asakura-Oosawa (AO) model Asakura and Oosawa (1954); Binder, Virnau, and Statt (2014) of colloid-polymer mixtures Poon (2002). In this model, colloidal particles and polymers are described as spheres of different sizes, with pairwise interactions. The ratio of particle sizes is $q\leq 1$. In the corresponding CG model, only the colloidal particles are included. If $q$ is sufficiently small then a CG model with pairwise interactions provides an exact description of the colloidal degrees of freedom Dijkstra, Brader, and Evans (1999). For larger $q$, only approximate CG models are available, because of effective interactions among the large particles that cannot be described by a pairwise additive potential Dijkstra, Brader, and Evans (1999); Lo Verso et al. (2006); Ashton and Wilding (2014). In this case – as in many other CG models which restrict to pairwise interactions Likos (2001); Ashton et al. (2011) – corrections to the CG model are required Ashton and Wilding (2014) in order to accurately estimate properties of the full (fine-grained) mixture.

For the specific case of the AO model, a method for calculation of these corrections was previously exploited by Dijkstra and co-workers Dijkstra and van Roij (2002); Dijkstra et al. (2006). They used a geometrical construction to analyse many-body interaction terms on the fly, during a CG simulation. Our method provides an alternative route to computation of these corrections: it provides an accurate characterisation of the critical point of the AO mixture. As expected Dijkstra and van Roij (2002); Dijkstra et al. (2006); Lo Verso et al. (2006), we find that many-body effective interactions for the colloids tend to reduce the tendency of the system to form clusters. We also analyse the nature and strength of the three-body and many-body effective interactions. In particular, for $q=\frac{2}{5}$, we show that a CG model with two- and three-body interactions can capture the critical point to high accuracy, where a two-body CG model fails.

The structure of the paper is as follows. Sec. II describes the general method, while Sec. III describes the AO model and details how the method is applied in this case. In Sec. IV we present numerical results. Then, Sec. V analyses the accuracy of the method, and its statistical uncertainties. Finally, Sec. VI summarises the conclusions and outlook. Additional technical detail is provided in the appendices.

II The Two-level Method

II.1 Outline

This Section presents our method in general terms, using a schematic notation. The specific application that we consider is a mixture of large and small particles, for which full definitions are given in Sec. III.

Consider a system with many degrees of freedom, and suppose that we are given a coarse-grained model which describes the same physical system in terms of a smaller set of (coarse) variables. In our application, the full system consists of a mixture of large and small particles, in the grand canonical ensemble. The coarse degrees of freedom are the number of large particles $N$, and their positions $(\bm{R}_{1},\bm{R}_{2},\dots,\bm{R}_{N})$. The other degrees of freedom are the number of small particles $n$, and their positions $(\bm{r}_{1},\bm{r}_{2},\dots,\bm{r}_{n})$. In our schematic notation, we denote a configuration of the coarse system as $\mathcal{C}=(N,\bm{R}_{1},\bm{R}_{2},\dots,\bm{R}_{N})$ and a configuration of the full system as $(\mathcal{C},\mathcal{F})$ with $\mathcal{F}=(n,\bm{r}_{1},\bm{r}_{2},\dots,\bm{r}_{n})$. This full system is the FG system, whose degrees of freedom comprise both $\mathcal{C}$ and $\mathcal{F}$. See Fig. 1(a).

Given some configuration of this FG system, the energy is $E(\mathcal{C},\mathcal{F})$ and the associated Boltzmann distribution is

[TABLE]

where $\beta=1/(k_{\rm B}T)$ is the inverse temperature (scaled by Boltzmann’s constant $k_{\rm B}$), and $Z$ is the partition function.

The aim of the method is to compute averages with respect to $p_{\rm f}$. We restrict to averages of observable quantities that depend only on the coarse degrees of freedom. The average of such a quantity $A$ is denoted by

[TABLE]

The definition of the integration measure $\mathrm{d}\mathcal{C}\mathrm{d}\mathcal{F}$ depends on the system of interest. The marginal probability density for the coarse variables is

[TABLE]

Hence $\langle A\rangle_{\rm f}=\int p_{\rm m}(\mathcal{C})A(\mathcal{C})\mathrm{d}\mathcal{C}$.

A configuration $\mathcal{C}$ of the coarse-grained system has energy $E_{\rm c}(\mathcal{C})$. The Boltzmann distribution for the coarse system is

[TABLE]

and the associated equilibrium average is

[TABLE]

For a given fine system, the identification of suitable coarse degrees of freedom $\mathcal{C}$ and the inference of the energy function $E_{\rm c}(\mathcal{C})$ are both part of the coarse-graining procedure. A related procedure is the generation of representative configurations of the full (FG) system, starting from coarse configurations $\mathcal{C}$; this is known as reverse-mapping, see for example Ref. Spyriouni et al., 2007. See Fig. 1(a).

A great deal of effort Likos (2001); Behler and Parrinello (2007); Bartók et al. (2010); Ouldridge, Louis, and Doye (2011); Mladek et al. (2013); Pak and Voth (2018) has gone into developing energy functions $E_{\rm c}$ such that coarse averages such as $\langle A\rangle_{\rm c}$ provide good estimates of the corresponding fine average $\langle A\rangle_{\rm f}$. This obviously requires that $p_{\rm c}\approx p_{\rm m}$. Our method addresses the following question: Suppose that $E$ and $E_{\rm c}$ are given, with $p_{\rm c}$ close (but not equal) to $p_{\rm m}$. How should one obtain an accurate estimate of $\langle A\rangle_{\rm f}$, at reasonable computational cost?

The TL method can achieve this, as shown schematically in Fig. 1(b). It combines computations on the CG model with a reverse-mapping procedure that is used to assess the accuracy of the this model. We first estimate $\langle A\rangle_{\rm c}$ by classical Monte Carlo (MC) sampling, which is assumed to have a relatively low computational cost. (In mathematics, this would be called a Markov chain Monte Carlo (MCMC) method.) This method is used to generate $\mathcal{N}_{\rm c}$ representative samples from (4) which should be identically distributed but need not be independent. These are denoted by $(\mathcal{C}_{1},\mathcal{C}_{2},\dots,\mathcal{C}_{\mathcal{N}_{\rm c}})$. Define the estimator

[TABLE]

This is a consistent estimator of $\langle A\rangle_{\rm c}$ but it is not an accurate estimate of the quantity of interest, which is $\langle A\rangle_{\rm f}$. To quantify the difference between these two averages we make a computational estimate of

[TABLE]

where $w(\mathcal{C})=p_{\rm m}(\mathcal{C})/p_{\rm c}(\mathcal{C})$ is a reweighting factor. This is achieved by drawing $\mathcal{N}_{\rm f}$ samples from those that were used in (1). With slight abuse of notation we denote these by $(\mathcal{C}_{1},\mathcal{C}_{2},\dots,\mathcal{C}_{\mathcal{N}_{\rm f}})$. Then one estimates (7) as

[TABLE]

where $\hat{w}(\mathcal{C}_{i})$ is an estimate of the ratio $p_{\rm m}(\mathcal{C}_{i})/p_{\rm c}(\mathcal{C}_{i})$ that is obtained by a free energy computation using Jarzynksi’s equality (see Section III.4, this step has significant computational cost but is tractable in the systems of interest). Finally, one has an estimator for the fine average:

[TABLE]

We show below that as $\mathcal{N}_{\rm c},\mathcal{N}_{\rm f}\to\infty$ then

[TABLE]

II.2 General properties

From (3), one identifies $(-1/\beta)\log p_{\rm m}(\mathcal{C})$ as the free energy of the $\mathcal{F}$ variables, for a given coarse configuration $\mathcal{C}$. Hence $E_{\rm c}(\mathcal{C})=(-1/\beta)\log[Z_{\rm c}p_{\rm c}(\mathcal{C})]$ is an estimate of this free energy, and the correction (9) is analogous to the free-energy perturbation formula of Zwanzig Zwanzig (1954). In general, this perturbation formula suffers from large errors, except when the correction is small. Our TL method provides a method for computing the weights $\hat{w}$ so as to obtain accurate estimates for $\langle A\rangle_{\rm f}$.

It is natural to compare the TL method with simple MC algorithms that operate directly on the fine system. In cases where coarse-graining is useful, one should expect a separation of length scales between coarse and fine degrees of freedom. Often, this means that direct MC calculations on the fine system are frustrated by the disparate length scales. For example, it is likely that MC moves should involve co-ordinates that change only on the fine scale, but exploration of configuration space requires changes over the (much larger) coarse scale. The result is that it takes many such moves to relax (or decorrelate) the configurations of the system. In mathematical parlance, the MCMC chain for the FG system is expected to mix very slowly.

The TL method avoids this problem: the central idea is illustrated in Fig. 1(b). Note that the estimator ${\hat{A}}_{\rm c}$ only involves coarse degrees of freedom, so it does not suffer from the separation of length scales. Computation of ${\hat{\Delta}}$ requires MC computations to be performed on the fine degrees of freedom but it does not require that the coarse degrees of freedom are relaxed during this process. Hence, this computation does not suffer from the (very slow) mixing in the FG system. In cases where the CG model is accurate, this can lead to a strong improvement in performance.

Given a FG system, the TL method can be applied using any CG model. The results are (in principle) correct in all cases. This feature – which is standard within methods based on histogram reweighting, importance sampling, and free-energy perturbations Bruce and Wilding (2003); Zwanzig (1954) – allows enormous flexibility in the choice of CG model. However, if the CG model is a poor representation of the fine one, the method is not efficient from a numerical point of view. In this case, the weights $\hat{w}_{i}$ acquire large fluctuations, which leads to large statistical errors on $\hat{\Delta}$ (that is, large error bars). This is an in-built diagnostic of the quality of the CG model, within the method. The computed weights can also be used to refine CG models, to improve their accuracy.

III Asakura-Oosawa model system

III.1 Fine-grained system

As an application of the TL method, we consider a model of a colloid-polymer mixture, as studied by Asakura and Oosawa Asakura and Oosawa (1954) and by Vrij Vrij (1976) It consists of large and small spherical particles in three spatial dimensions, with diameters $\sigma_{\rm B}$ and $\sigma_{\rm S}$ respectively. The behaviour of this model is reviewed in Ref. Binder, Virnau, and Statt, 2014, see also Fig. 2. Its phase behaviour has been studied extensively, both computationally Dijkstra, Brader, and Evans (1999); Dijkstra et al. (2006) and theoretically Gast, Hall, and Russel (1983); Lekkerkerker et al. (1992); Schmidt et al. (2002); Brader, Evans, and Schmidt (2003); Oversteegen and Roth (2005). The large particles are hard spheres that represent the colloid. The small particles are schematic representations of (ideal) polymer chains and their size is given by the polymer radius of gyration. These spheres are free to overlap with each other, because the excluded volume of a polymer is very small in comparison with its radius of gyration. In colloid-polymer mixtures, the effect of the polymers is to promote clustering of the colloidal particles, which leads eventually to phase separation and/or crystallisation Poon (2002). The AO model captures both these features Binder, Virnau, and Statt (2014). We focus here on fluid phases.

The particles occupy a periodic box in three dimensions, of size $L\times L\times L$ . The large particles are hard: the potential energy $U_{\rm f}(\mathcal{C})=\infty$ if two large particles overlap or if a large particle overlaps with a small one; otherwise $U_{\rm f}(\mathcal{C})=0$. This means in particular that small particles are free to overlap with each other. The chemical potentials of the large and small particles are $\mu_{\rm B},\mu_{\rm S}$. The Boltzmann distribution for this system is similar to (1):

[TABLE]

where $\Xi_{\rm f}$ is the grand partition function. Computation of averages like (2) requires a definition of the integration measure $\mathrm{d}\mathcal{C}\mathrm{d}\mathcal{F}$: this is chosen such that the average in (2) is equal to an average over the grand canonical ensemble. It includes sums over $N,n$ and also a combinatorial factor of $N!n!$ that accounts for particle indistinguishability. Details are given in Appendix A.1.

In place of the parameter $\beta\mu_{\rm S}$, it is convenient to work with the volume fraction of the small particles $\eta_{\rm S}^{\rm r}$ in a (nominal) reservoir at chemical potential $\mu_{\rm S}$, that is $\eta_{\rm S}^{\rm r}=\pi\mathrm{e}^{\beta\mu_{\rm S}}/6$. The size ratio between large and small particles is

[TABLE]

Since the small particles do not interact with one another, the integral in (3) can be computed, by summing over the number of small particles and integrating their positions. Following Ref. Dijkstra, Brader, and Evans, 1999, we show in Appendix A.1 that for configurations $\mathcal{C}$ without overlapping particles the result is

[TABLE]

where ${\cal V}_{\rm a}(\mathcal{C})$ is the total volume that is accessible to the (centres of) the small particles.

We summarise the dimensionless model parameters, which are the colloid size ratio $q$, the large-particle chemical potential $\beta\mu_{\rm B}$, the volume fraction of the small particles (in the reservoir) $\eta_{\rm S}^{\rm r}$, and the size of the simulation box $L/\sigma_{\rm B}$.

III.2 Coarse-grained model

The coarse system only has large particles. These still have (repulsive) hard cores but they also interact by a pairwise additive interaction with the Asakura-Oosawa (AO) pair potential. That is, for configurations $\mathcal{C}$ without overlapping particles one has Binder, Virnau, and Statt (2014)

[TABLE]

with

[TABLE]

If there are overlapping particles then $U_{\rm c}=\infty$. The corresponding Boltzmann distribution is

[TABLE]

where $\Xi_{\rm c}$ is the coarse partition function. Comparing (13,16), one sees that the coarse model corresponds to a pairwise additive approximation for the accessible volume ${\cal V}_{\rm a}$. The model is parameterised by the same dimensionless parameters as the fine system ($q,\beta\mu_{\rm B},\eta_{\rm S}^{\rm r},L/\sigma_{\rm B}$).

The behaviour of this system depends on the size ratio $q$. In particular, for $q\leq(2/\sqrt{3})-1\approx 0.1547$, a geometrical construction can be used to show that the pairwise-additive approximation for ${\cal V}_{\rm a}$ is exact and $p_{\rm c}(\mathcal{C})=p_{\rm m}(\mathcal{C})$ for all coarse configurations $\mathcal{C}$ Binder, Virnau, and Statt (2014). For larger values of $q$, the effective interaction in the fine system includes 3-body and higher contributions, which means that

[TABLE]

where the sums run over all the large particles, $\Delta\mu$ is a constant, and $V_{3},V_{4},\dots$ are many-body interaction potentials. For the specific case of the AO model, these many-body terms can be computed (numerically) using a geometrical construction, as shown in Refs. Dijkstra and van Roij, 2002; Dijkstra et al., 2006. For the purposes of this work we note that $p_{\rm m}\neq p_{\rm c}$ so computations on the coarse system (16) do not provide accurate estimates of properties of the fine system. Hence the method of Sec. II is applicable.

The estimators (6,8) require that we generate representative samples from (16). This is achieved by grand canonical MC simulations of the coarse system Frenkel and Smit (2001).

III.3 Observable quantities

The phase diagram for the AO system is sketched in Fig. 2(a), as a function of $\eta_{\rm S}^{\rm r}$ and the large particle number density $\rho_{\rm B}$, for a representative value of $q$. For small $\eta_{\rm S}^{\rm r}$ there is a fluid phase at low to moderate density, and a crystalline phase at high density. On increasing $\eta_{\rm S}^{\rm r}$ at moderate $\rho_{\rm B}$, the crystalline phase becomes stable, but the fluid is still accessible as a metastable phase. On further increasing $\eta_{\rm S}^{\rm r}$ within the metastable fluid, there is a critical point at $(\eta_{\rm S}^{\rm r})_{\rm c}$, and for $\eta_{\rm S}^{\rm r}>(\eta_{\rm S}^{\rm r})_{\rm c}$ one observes coexistence between phases with many large particles (“colloidal liquid”) and few large particles (“colloidal vapour”). This work concentrates on the behaviour near the critical point. In this regime the system depends very sensitively on interaction parameters, so it provides a challenging test case for the TL method. In our grand-canonical MC simulations of the coarse system, the fluid is sufficiently metastable that we do not observe crystallisation (see also Appendix B).

In the following we focus on the critical point. As a quantity of interest, we consider the histogram of the number of large particles. That is, we compute the probability to have exactly $N_{\rm b}$ large particles, which is

[TABLE]

This is of the required form (2). The distribution $P_{\rm b}$ has a characteristic form at the critical point Wilding (1995). In the vicinity of this point, it depends strongly on the model parameters. Hence the accurate determination of $P_{\rm b}$ is a stringent test of our method and allows us to locate the critical point of the fine system.

A previous study Lo Verso et al. (2006) has characterised this critical point by direct Monte Carlo computations on the fine system. However, despite the use of advanced MC methods Vink and Horbach (2004), that work was limited to relatively large $q\geq\frac{2}{5}$, otherwise the scale separation between the large and small particles makes such computations extremely expensive. Fig. 2 shows a typical configuration of the FG system near the critical point, for a small size ratio $q=\frac{2}{13}$. This picture illustrates that particle insertion MC moves for the large particles would be extremely challenging in this regime. Here we are able to obtain accurate results for the fine system at $q=\frac{1}{4}$. This allows us to estimate effects of three-body interactions, which are not captured by the coarse model, although they could be analysed following Ref. Dijkstra et al., 2006. Our method is also accurate at $q=\frac{2}{13}\approx 0.15$: this case is useful as a test of the method but the coarse model is already exact at that state point so it does not yield any new physical insight.

III.4 Jarzynski integration and estimation of $P_{\rm b}$

Computation of the weights $\hat{w}$ in (8) involves grand canonical MC simulations where the large particles are held fixed while small particles are inserted and removed. The system is divided into cells of size $\sigma_{\rm S}$ so that there are $n_{\rm cell}$ cells in total: this facilitates the computation of whether an attempted small-particle insertion is successful. We define a Monte Carlo sweep (MCS) as a sequence of $n_{\rm cell}$ attempted MC updates, where each update is chosen to be either a particle insertion or a particle removal, with equal probability.

The weight factors $\hat{w}(\mathcal{C}_{i})$ are computed in two stages. We first calculate non-normalised weights $\hat{W}(\mathcal{C}_{i})$ by a Monte Carlo (Jarzynski) integration method Hummer and Szabo (2001); Neal (2001). The procedure is described below where we also explain that for any fixed $\mathcal{C}$ then

[TABLE]

The average in this equation is an expectation value with respect to the Monte Carlo (Jarzynski) integration (see below). Then we take

[TABLE]

Observe from (13,16) that (19) is equal to $(\Xi_{\rm f}/\Xi_{\rm c})\cdot(p_{\rm m}(\mathcal{C})/p_{\rm c}(\mathcal{C}))$. Taking the average of (20) and comparing with (13,16), one has for large $\mathcal{N}_{\rm f}$ that

[TABLE]

where we used that the denominator of (20) converges for large $\mathcal{N}_{\rm f}$ to the expectation of (19) with respect to $p_{\rm c}$, which is $(\Xi_{\rm f}/\Xi_{\rm c})$. [The approximate equality in (21) is accurate up to a small correction associated with the difference between the denominator of (20) and its average, see Sec. V, below.]

We now describe the computation of the weights $\hat{W}$. Given some $\mathcal{C}$, we identify the grand partition function for the small particles as

[TABLE]

To compute $\hat{W}(\mathcal{C})$, we hold the large particles fixed in position and we perform grand-canonical MC for the small particles. We first choose some chemical potential $\mu_{0}$ which is small enough that the partition function $\Xi_{\rm S}(\mu_{0},\mathcal{C})$ can be estimated directly from a grand canonical simulation of the small particles: it is equal to the reciprocal of the probability that $n=0$.

Next we perform an MC simulation with $M$ MCS during which the small particle chemical potential increases from $\mu_{0}$ to $\mu_{\rm S}$. Let the small-particle chemical potential after $t$ MCS be $\mu(t)$, and the number of small particles at this time be $n(t)$. Then define

[TABLE]

The parameter $M$ and the dependence of $\mu$ on $t$ are arbitrary at this stage and can be chosen to optimise the accuracy of the method. In our simulations, $\mu$ increases stepwise. Let $\Delta\mu_{j}$ be the change in $\mu$ on the $j$th step and $n_{j}$ be the number of small particles in the system when that step takes place. Then

[TABLE]

where $K$ is the number of steps. (A simple choice is to make one step per MCS so that $M=K$, but other choices are possible.) A consistent choice of the weight is then $\hat{W}(\mathcal{C})=\Xi_{\rm S}(\mu_{0},\mathcal{C})\mathrm{e}^{\beta{\cal I}(\mathcal{C})+\beta U_{\rm c}(\mathcal{C})}$. Then Crooks’ fluctuation formula Crooks (2000) (analogous to Jarzynski’s equality Jarzynski (1997)) can be used to show that (19) holds. A short derivation is given in Appendix A.2, see also Refs. Hummer and Szabo, 2001; Neal, 2001; Oberhofer and Dellago, 2009. Alternatively we may repeat the same computation $m$ times (always with the same $\mathcal{C}$) and take

[TABLE]

where ${\cal I}_{p}(\mathcal{C})$ is the value of (24) obtained on the $p$th iteration of the computation. Finally, using (25) in (20) gives the weight factors $\hat{w}_{i}$.

We summarise the parameters of the TL method. The numbers of configurations used for the coarse and fine estimates are $\mathcal{N}_{\rm c},\mathcal{N}_{\rm f}$. These should both be large in order to have accurate results. In addition, the Jaryznski integration depends on $M$ which is the total number of MCS in a Jarzysnki run, and $m$ which is the number of Jarzynski runs for each $\mathcal{C}$. The computation of $\hat{\Delta}$ is the dominant computational cost in the method, which scales as $mM\mathcal{N}_{\rm f}$. The choice of these parameters is discussed in Sec. V. The other choice that must be specified is the annealing schedule: that is, the dependence of $\mu$ on $t$ in (23), as well as the initial chemical potential $\mu_{0}$ and the number of steps $K$. The choice of this schedule is discussed in Appendix B.

IV Numerical results

IV.1 Critical behaviour

As anticipated in Sec. III.3, we focus on the probability distribution of $N$, the number of large particles in the system. We analyse this distribution based on the theory of liquid-vapour critical points, but neglecting field-mixing effects Wilding (1995). At this level of theory, the distribution of $N$ has a characteristic (universal) form at the critical point, with two peaks of equal height, separated by a trough. The ratio of probability between the peak and trough is 0.46, in contrast to the phase-coexistence regime for which this ratio scales with system size as $\mathrm{e}^{-\kappa L^{2}}$, where $\kappa$ is proportional to the surface tension. The two peaks correspond to particle densities whose difference scales as $L^{-\beta/\nu}$ where $\beta,\nu$ are the critical exponents for the 3d Ising universality class, so that $\beta/\nu\approx 0.52$, see for example Ref. Wilding, 1995. Hence the variance of the distribution decreases as $L^{-2\beta/\nu}$, as expected at liquid-vapour critical point. (These critical fluctuations are much larger than a comparable three-dimensional system with short-ranged correlations, where the variance decreases as $L^{-3}$).

In order to characterise the shape of the distribution, we define a scaling variable $X$ associated with the number of large particles by subtracting the mean of $N$ and dividing by the standard deviation:

[TABLE]

The probability density for this variable is denoted by $P_{L}(X)$, where $L$ indicates the system size.

IV.2 Size ratio $q=\frac{1}{4}$

For $q=(1/4)$, direct simulation of the FG system is challenging, especially in the grand-canonical ensemble. However, the TL method is applicable. Fig. 3 shows results obtained using this method, in a system of size $L=5\sigma_{\rm B}$. Fig. 3(a) shows that the critical point in this system occurs at $\eta_{\rm S}^{\rm r}\approx 0.4080$. Taking this value and making an appropriate choice of $\mu_{\rm B}$, the distribution of $N$ matches well to the universal critical distribution.

Fig. 3(b) shows results at $\eta_{\rm S}^{\rm r}\approx 0.4025$, for both coarse and fine models. One sees that the coarse model is very close to its critical point, but the corresponding distribution for the fine model does not match the critical form. In the fine system, the bimodality of the distribution is less pronounced than the critical form, indicating that this system is in the one-phase regime, for this value of $\eta_{\rm S}^{\rm r}$. The physical interpretation of this result is that three-body (and higher) interactions in the fine system act to reduce the attractive forces, which pushes the critical point to a larger value of $\eta_{\rm S}^{\rm r}$. This shift is $\Delta\eta_{\rm S}^{\rm r}\approx 0.006$ for these parameters. This is a small effect but our method is accurate enough to resolve it. The qualitative difference between the FG distributions in Figs. 3(a,b) emphasises that critical properties are very sensitive to model parameters: the values of $\eta_{\rm S}^{\rm r}$ differ by less than $2\%$.

Note also that estimating the critical point by this method is only accurate up to an error that decays as a power law in $L$. This error has contributions from field-mixing effects and from corrections to scaling Wilding (1995). For the systems considered here, this correction might be comparable with $\Delta\eta_{\rm S}^{\rm r}$ itself. However, we expect a similar correction for both coarse and fine systems, in which case our measured $\Delta\eta_{\rm S}^{\rm r}$ should be a reasonable estimate of the difference in the critical points between coarse and fine systems.

For the system shown in Fig. 3(a), a histogram of the reweighting factors $\hat{w}_{i}$ is shown in Fig. 3(c). If the coarse and fine models matched exactly and the Jarzynski integration was perfect then all weight factors would be equal to unity. We see that some weights are significantly larger than $1$: these are predominately coming from configurations whose probabilities are underestimated by the coarse model. Such configurations must be given larger weight when estimating observable quantities for the fine model, and this is the effect of the correction $\hat{\Delta}$.

IV.3 Size ratio $q=\frac{2}{5}$

We now turn to a size ratio $q=(2/5)$. This system was analysed by direct simulation in Ref. Lo Verso et al., 2006, using the grand-canonical ensemble, with cluster moves Vink and Horbach (2004) that remove small particles at the same time as inserting large ones (and vice versa). As in this work, the critical point was estimated by comparison with the universal order-parameter distribution: they found $\eta_{\rm S}^{\rm r}=0.5215\pm 0.0001$, in a system of size $L=10\sigma_{\rm B}$. We note that the TL method does not outperform these previous methods for this moderate value of $q$: the TL method is designed to be effective for smaller $q$, which is the most challenging regime. For example, the next section considers $q=(2/13)$.

Ref. Lo Verso et al., 2006 also shows that three-body interactions are significant for the size ratio $q=(2/5)$, in that the critical points differ significantly between the fine model and the coarse model with energy (14). We show in Appendix C that the differences between the coarse and fine systems are too large for the TL method to be able to correct the coarse-graining error in that case We also explain how this problem can be avoided. Specifically, we define a new coarse model whose energy is

[TABLE]

with $U_{\rm c}$ given by (14) and

[TABLE]

where $R_{ab}=|\bm{R}_{a}-\bm{R}_{b}|$ is the distance between particles $a$ and $b$. That is, $V_{3}$ is a 3-body interaction potential which we have chosen to represent as a function of the three relevant interparticle distances [see Appendix C and recall also (17)].

We refer to the model with energy (27) as the (2+3)-body coarse model. The three-body interaction is repulsive: it weakens the effective attractive attraction between particles so the critical point for the (2+3)-body coarse model is at a higher volume fraction $\eta_{\rm S}^{\rm r}$ than the critical point of the two-body coarse model (14).

Results using the (2+3)-body coarse model within the TL method are shown in Fig. 4. In particular, Fig. 4(a) shows histograms of the large-particle density for coarse and fine models at $\eta_{\rm S}^{\rm r}=0.5174$, compared with the critical form. Two results are notable: first, the fine model matches well the critical distribution, so we conclude that this system is close to its critical point. Second, the fine model matches extremely accurately the coarse model. This indicates that the (2+3)-body coarse model is an almost-quantitative description of the fine system. This is confirmed by Fig. 4(b) which shows the distribution of the $\hat{w}_{i}$ which is peaked near $\hat{w}_{i}=1$, consistent with the fact that the coarse-grained model is very accurate. Compared with Ref. Lo Verso et al., 2006, the difference (less than $1\%$) in the estimated critical value of $\eta_{\rm S}^{\rm r}$ is presumably attributable to the smaller system size $L=6.4\sigma_{\rm B}$ used in this work.

Since the (2+3)-body coarse model gives an almost-exact representation of the critical point, we conclude that four-body and higher interactions are negligible at this state point. This is perhaps not surprising since the (average) volume fraction of the large particles is around $0.22$ at the critical point, and the effects of four-body and higher interactions are suppressed as the fourth power of the density. However, the result is promising in that it shows how low-order corrections to two-body interaction potentials can already give a significant improvement in accuracy.

IV.4 Size ratio $q=\frac{2}{13}$

As a third example system, we consider $q=(2/13)\approx 0.1538$ for which the two-body coarse model is exact Binder, Virnau, and Statt (2014). Results are shown in Fig. 5. They show that the FG and CG models do indeed have the same order-parameter distribution, as they must. There are many small particles (recall Fig. 2, which has over 20,000 of them) so direct simulation of this system would be extremely challenging.

Within the TL method it is known that $\langle\hat{\Delta}\rangle=0$ because the coarse model is exact. However, it is not easy to confirm this result numerically because the Jarzynski integration requires significant computational cost. Nevertheless, Fig. 5(b) shows that the distribution of weights $\hat{w}_{i}$ is peaked near unity, which means that the method is effective, even for this challenging test case. Fig. 5(a) confirms that the results are accurate. The width of the distribution of weights comes entirely from the Jarzynski integration. It is suppressed by choosing a large value of $M$, corresponding to a significant computational effort for the computation of each individual weight. The small error on $P(N)$ is ensured by combining a sufficiently narrow distribution of weights (large $M$), with an average over many of these weights (sufficiently large $\mathcal{N}_{\rm f}$). For details of the simulations see Appendix B. The scaling of the numerical error with respect to the algorithmic parameters is the subject of the next section.

The physical insight of Fig. 5 is that the TL method is applicable even when the number of small particles is very large, recall Fig. 2(b).

V Error estimation and choice of algorithmic parameters

V.1 Overview

This section analyses the convergence shown in (10). We will show that as $\mathcal{N}_{\rm c},\mathcal{N}_{\rm f}\to\infty$ then

[TABLE]

where the outer average (without subscript) is over all random aspects of the TL algorithm, which includes the generation of the original $\mathcal{N}_{\rm c}$ coarse samples, the selection of the $\mathcal{N}_{\rm f}$ samples to be used in (8), and the Jarzynski integration. The mean square error of the method is

[TABLE]

We will analyse the scaling of the error $h^{2}$ in the limit $\mathcal{N}_{\rm c},\mathcal{N}_{\rm f}\to\infty$.

Our analysis requires several assumptions, which are all satisfied in the physical situations of interest. Note first that the weighting factor $w(\mathcal{C})=p_{\rm m}(\mathcal{C})/p_{\rm c}(\mathcal{C})$ may be very large if $p_{\rm c}$ is very small. This happens if there are configurations of the FG model that are very rare in the CG model. We assume that the variance of $w(\mathcal{C}_{j})$ is finite, which requires that the CG model is not too inaccurate. In fact, it is useful to quantify the accuracy of the CG model in terms of the $\chi^{2}$-divergence Chen (2005); Agapiou et al. (2017) between $p_{\rm c}$ and $p_{\rm m}$, which is defined as

[TABLE]

This quantity is non-negative in general and equals zero only if the coarse model is perfect, that is $p_{\rm c}=p_{\rm m}$. Smaller values of $D_{\chi^{2}}$ correspond to coarse models that are more accurate. Our assumption is that

[TABLE]

We also require an assumption on the variance of the weight factors $\hat{W}$ generated by the Jarzynski integration method. As a minimal assumption, we suppose that this variance is finite, which means that

[TABLE]

The expected behaviour of the weights $\hat{W}$ is discussed below, but we expect that (33) holds in all situations of physical interest.

In addition to (32,33), we make several other assumptions, to simplify the analysis. In particular, we assume that all the samples (both coarse and fine) are independent of each other. In practical situations this assumption will not hold exactly, but we argue that our results still have the correct scaling with respect to the algorithmic parameters. For simplicity, we also assume throughout that $A$ is a bounded quantity, so that

[TABLE]

for some $a$, independent of $\mathcal{C}$.

V.2 Convergence of the mean

To establish (29), we first note from (6) that

[TABLE]

which follows by linearity of the average, because every $\mathcal{C}_{i}$ is distributed as $p_{\rm c}$. Similarly,

[TABLE]

The next step is to quantify the approximation in (21). To this end define $\xi=\Xi_{\rm f}/\Xi_{\rm c}=\langle\hat{W}(\mathcal{C}_{j})\rangle$ and also

[TABLE]

Clearly $\langle\epsilon\rangle=0$, and $\epsilon$ is an average of independent mean-zero variables. Also $\hat{W}(\mathcal{C}_{j})$ has a finite variance by (33), so that $\epsilon$ obeys a central limit theorem, and its typical value is therefore small as $\mathcal{N}_{\rm f}\to\infty$. Also from (20)

[TABLE]

which is a more precise statement of (21). Observe that $\sum_{j}\hat{w}(\mathcal{C}_{j})=1$ so (8) is identical to

[TABLE]

Hence, using again that the $\mathcal{C}_{j}$ are distributed as $p_{\rm c}$, one has

[TABLE]

with

[TABLE]

Using (38) and that $\epsilon$ is small, we can establish a bound on the average of $\hat{\mathcal{E}}$. Details are given in Appendix D. The result is

[TABLE]

With (9,35), this establishes (29).

We have emphasised throughout that the Jarzynski method for computing weights is useful because it minimises systematic errors. To be specific, the simplest alternative to this approach would be to use a standard thermodynamic integration to estimate $p_{\rm m}/p_{\rm c}$, with the same computational effort $M$ as was used in the Jarzynski method. The disadvantage of this approach is that one expects a systematic error on $\hat{\Delta}$ from the finite rate of integration, which would appear in (42) as a correction at ${\cal O}(M^{-1})$. (Recall that $M$ is the total number of MCS used in the Jarzynski integration.) In this case convergence of the method requires a joint limit $\mathcal{N}_{\rm c},\mathcal{N}_{\rm f},M\to\infty$, instead of the (simpler) limit $\mathcal{N}_{\rm c},\mathcal{N}_{\rm f}\to\infty$ required here.

V.3 Mean square error

This section estimates the the error defined in (30), and how it scales with the parameters $\mathcal{N}_{\rm c},\mathcal{N}_{\rm f},M$. Let

[TABLE]

Assuming that the $\mathcal{N}_{\rm c}$ samples used in $\hat{A}_{\rm c}$ are independent, it is immediate that

[TABLE]

Using (34), the numerator is bounded above by $a^{2}$, so $h_{\rm c}^{2}\to 0$ as $\mathcal{N}_{\rm c}\to\infty$. In practice, (5) uses a sequence of samples from a MC simulation, so that samples which are nearby in the sequence are not independent. However, these correlations will be short-ranged, as long as the MC simulation is converged. In this case $h_{\rm c}^{2}$ is larger than (44) by a constant factor but the scaling with $\mathcal{N}_{\rm c}$ is the same.

For the $\mathcal{N}_{\rm f}$ samples used in (8), it is again convenient to assume that these samples are independent of each other and of the coarse samples. The mean square error of $\hat{\Delta}$ is

[TABLE]

Writing $\hat{\Delta}=(\hat{\Delta}-\langle\hat{\Delta}\rangle)+\langle\hat{\Delta}\rangle$ and using (42) yields

[TABLE]

We show below that $h_{\rm f}^{2}=O(\mathcal{N}_{\rm f}^{-1})$ so it is consistent to drop the correction term. Using (39) we then obtain the leading order result

[TABLE]

Substituting for $\hat{w}$ with (38), the leading-order behaviour can be obtained by dropping all terms proportional to $\epsilon$. Hence

[TABLE]

It follows from (34) that for any random variable $B$ then $\mathrm{Var}[(A(\mathcal{C}_{j})-\langle A\rangle_{\rm f})B(\mathcal{C}_{j})]\leq a^{2}\langle B(\mathcal{C}_{j})^{2}\rangle$. Hence (48) reduces to

[TABLE]

The physical interpretation of (49) is that the scaling of the error is controlled by the variance of the weights $\hat{W}$. There are two contributing factors to this variance, which are the randomness in the Jarzynski integration and the fact that $\mathcal{C}_{i}$ is distributed as $p_{\rm c}$.

We split the error estimate into two pieces, one from each source of randomness. Assume first that the reweighting factors $\hat{w}$ are obtained exactly so $\hat{w}_{i}=p_{\rm m}(\mathcal{C}_{i})/p_{\rm c}(\mathcal{C}_{i})$. In this case all randomness comes from the $\mathcal{C}_{i}$. Hence (49) reduces to

[TABLE]

The case where $p_{\rm c}=p_{\rm m}$ is instructive, it corresponds to a perfect coarse model so that $D_{\chi^{2}}=0$ and $h_{\rm f,0}^{2}=0$. Note however that the inequality in (50) is valid only if the weights $\hat{w}$ are exact.

To account for the randomness in these weights, recall that the average of $\hat{W}(\mathcal{C})$ obeys (19), but (48) also depends on the fluctuations of this quantity. To estimate the variance of $\hat{W}(\mathcal{C})$, note that ${\cal I}$ in (24) is a sum of weakly correlated increments, so one expects its statistics to be ruled by a central limit theorem Neal (2001). For this reason, we assume that for any fixed $\mathcal{C}$, the sum ${\cal I}$ has a Gaussian distribution with variance $\Lambda/M$. We assume that $\Lambda$ is independent of $\mathcal{C}$, for convenience.

Recall from (25) that we obtain the weight $\hat{W}$ by averaging over $m$ independent realisations of the Jarzynski integration. If $m=1$ one sees that $\log\hat{W}(\mathcal{C})$ differs from ${\cal I}$ by an additive constant and is therefore Gaussian with variance $\Lambda/M$. It follows that $\langle\hat{W}(\mathcal{C})^{2}\rangle_{\rm J}=\langle{W}(\mathcal{C})\rangle_{\rm J}^{2}\exp(\Lambda/M)$. For $m>1$ the variance of $\hat{W}$ is suppressed by a factor of $m$, which leads to

[TABLE]

Hence

[TABLE]

Using this result in (49) yields

[TABLE]

with

[TABLE]

The key point here is that even if the coarse model is perfect ($D_{\chi^{2}}=0$), there is randomness in the Jarzynski integration which generates statistical errors. From (50,53,54) one sees that indeed $h_{\rm f}^{2}=O(\mathcal{N}_{\rm f}^{-1})$, as asserted above.

Combining all these ingredients and recalling our various assumptions, we find the scaling of the mean square error,

[TABLE]

with

[TABLE]

V.4 Implications for optimisation of the TL method

The results (55,56) establish useful properties of the TL algorithm. The computational effort for the Jarzynski integration is proportional to the product $mM\mathcal{N}_{\rm f}$. We infer that the optimal parameter choice is to take $\Lambda/M=O(1)$ in order that the exponential factor in (56) is not too large, and then to devote the remaining computational effort towards larger $\mathcal{N}_{\rm f}$. The benefits of increasing $m$ are less than those of increasing $\mathcal{N}_{\rm f}$ or $M$, so the choice $m=1$ is likely optimal. However, it may be convenient in practice to take some other small integer since (for example) increasing $m$ can be achieved by additional parallel processing while increasing $M$ requires longer simulation (wall) time. The parameters of the Jarzynski integration (particularly the annealing schedule) should be chosen to minimise $\Lambda$ (at fixed $M$). It is obviously desirable that $D_{\chi^{2}}$ be small – this simply means that the coarse model should be as accurate as possible. In fact, $D_{\chi^{2}}$ can be related to the exponential of a free energy difference between $p_{\rm m}$ and $p_{\rm c}$, via (31). Hence we expect that this quantity will scale exponentially in the system size – very accurate coarse models are required if the method is to be effective in very large systems.

The factors of $D_{\chi^{2}}$ and $\mathrm{e}^{\Lambda/M}$ in (56) both indicate that the method can suffer from exponentially growing errors. However, the numerical results of Sec. IV show that the method is still practical in realistic soft-matter systems. The fact that $D_{\chi^{2}}$ and $\Lambda$ may be large in practice indicate that the majority of the computational effort in this method should be spent on the computation of $\hat{\Delta}$, which requires that both $M$ and $\mathcal{N}_{\rm f}$ should be large. By contrast, we found that the error $h_{\rm c}$ coming from $\mathcal{N}_{\rm c}$ tends to be small, in the systems analysed here.

VI Conclusion

We have presented a two-level method for correction of coarse-graining errors and we have demonstrated it using the AO mixture. We have found that three-body interactions shift the critical point to higher $\eta_{\rm S}^{\rm r}$, for $q=(2/5)$ and $q=(1/4)$. For $q=(2/5)$ we characterised the three-body interactions, which have a significant effect on the critical point, and also we showed that four-body and higher interactions are negligible, at criticality.

We also analysed the errors associated with the method, in Sec. V. The first important result of that section is that the method is accurate as $\mathcal{N}_{\rm c},\mathcal{N}_{\rm f}\to\infty$ but does not require any limit of large $M$. That is, we do not require that the individual weight factors $\hat{w}_{i}$ are accurate measures of $p_{\rm m}/p_{\rm c}$. Instead we require only that these weights have the right average, which is ensured by (21). The results (55,56) show how the parameters of the method should be chosen, as discussed in Sec. V.4.

As noted above, the general two-level method is similar in spirit to free-energy perturbation theory and importance sampling. As such, it has a range of potential applications. Its main strength is that direct simulation of the whole system is never required. For example, in the system presented here, the large particles remain fixed while the small particles are introduced, which avoids the difficulty of moving (or inserting) large particles in a system like Fig. 2. If one attempted direct MC simulation of the full system then the error (55) would be inversely proportional to the total number of MCS, but it would be proportional to the number of MCS required to generate an independent sample of the full system. This number is expected to be extremely large in practice. (In particular, it would likely be much larger than $\Lambda$.) On the other hand, direct MC simulation does not suffer from problems associated with large values of $D_{\chi^{2}}$, which can hamper the TL method. That is, the TL method requires rather accurate coarse descriptions in order to be successful (see also Fig. 6 in Appendix C). Nevertheless, we have shown that information from Jarzynski integration can be used to learn about corrections to coarse models, including three-body interactions. Combination of the TL method with machine-learning methods Behler and Parrinello (2007); Bartók et al. (2010); Cheng et al. (2019) seems to be a promising future direction.

Of course, there are many strategies for direct MC simulation. In particular there are tailored algorithms for mixtures of large and small particles, including the cluster moves used in Ref. Lo Verso et al., 2006, other cluster moves Liu, Wilding, and Luijten (2006), and staged (grand-canonical) particle insertion Ashton and Wilding (2011). These methods reduce the number of MCS required to generate an independent sample. For the specific case of the AO model, the method of Ref. Dijkstra et al., 2006 is also available. A direct comparison with these methods is tricky because the computational efficiency depends strongly on details of the implementation. However, the results shown here are promising, especially given the relatively high volume fractions of the small particles. We also emphasise that the TL method is generic, while these other algorithms have been tailored to the specific problems of interest.

Another strength of the TL method is that while the computational effort of the Jarzynski integration is proportional to $mM\mathcal{N}_{\rm f}$, this effort is easily split into $m\mathcal{N}_{\rm f}$ independent simulations, each involving $M$ MCS. This means that the method is well-suited to parallel computation, in contrast to direct MC simulation where very long individual simulations may be required.

In the future, we intend to test the TL method on mixtures of large and small hard spheres, for which accurate coarse models can be computed theoretically Dijkstra, van Roij, and Evans (1998, 1999); Ashton et al. (2011). The method will allow analysis of phase behaviour in these systems, and also quantification of the effects of many-body interactions Ashton and Wilding (2014). Similar questions arises in systems where the large particles are not spherical, in which case the depletion interaction has a more complicated character Odriozola, Jimenez-Angeles, and Lozada-Cassou (2008); van Anders et al. (2014); Law et al. (2016) and can lead to new phase behaviour van Anders et al. (2014); Ashton, Jack, and Wilding (2015). There are also potential applications of the TL method to other soft matter systems including star polymers and dendrimers, for which accurate coarse-grained models are available Jusufi, Watzlawek, and Löwen (1999); Likos (2001); Marzi, Likos, and Capone (2012), and whose phase behaviour is rich and complex Lenz et al. (2012); Wilding and Sollich (2014).

Acknowledgements.

This work was supported by the Leverhulme Trust through research project grant number RPG–2017–203.

Appendix A Additional theoretical details

A.1 Grand canonical ensemble

To analyse these fluids in the grand canonical ensemble requires a definition of the integration measure in (2), which includes a sum over the numbers of particles and an integral over their positions. We take

[TABLE]

where all particle positions are integrated over the box $[0,L]^{3}$. This formula is consistent with (11), recall that $\mathcal{C}=(N,\bm{R}_{1},\dots,\bm{R}_{N})$ and $\mathcal{F}=(n,\bm{r}_{1},\dots,\bm{r}_{n})$ and note that the zero of chemical potential has been chosen so that the thermal de Broglie wavelength does not appear, as usual in simulation studies. Note also that the definition of the integration measure $\mathrm{d}\mathcal{C}$ includes the factor of $\frac{1}{N!}$ that accounts for indistinguishability of particles, and similarly for $\mathrm{d}\mathcal{F}$.

Also, normalisation of $p_{\rm f}$ requires

[TABLE]

Similar formulae hold for $\langle A\rangle_{\rm c}$ and $\Xi_{\rm c}$.

To compute $p_{\rm m}(\mathcal{C})$ from (3,11) one notes that the integrations over the particle positions can be performed independently and that each of them yields a factor of ${\cal V}_{\rm a}$, which is the volume accessible to a single small particle. For configurations $\mathcal{C}$ where no large particles are overlapping one finds

[TABLE]

Recognising the sum as the Taylor expansion of an exponential and using the definition of $\eta_{\rm S}^{\rm r}$ yields (13).

A.2 Jarzynski average (19)

We outline the derivation of (19), following Ref. Crooks, 2000. As in (24), the chemical potential increases in a sequence of $K$ steps. Hence we consider a sequence of $K+1$ values of the chemical potential, $(\mu_{0},\mu_{1},\dots,\mu_{K})$, with $\mu_{K}=\mu_{\rm S}$, the value of the chemical potential for the fine system. Also $\Delta\mu_{k}=\mu_{k}-\mu_{k-1}$. Between steps $k$ and $k+1$ one inserts and removes small particles using an MC algorithm that obeys detailed balance with respect to a probability distribution

[TABLE]

which is similar to (11) with $\mu_{\rm S}$ replaced by $\mu_{k}$. To establish (19) from (25) it is sufficient to show that

[TABLE]

where the average is with respect to this random MC algorithm and the small-particle partition function of (22) can be alternatively expressed as

[TABLE]

Let the sequence of small-particle configurations when the steps take place be $(\mathcal{F}_{1},\mathcal{F}_{2},\dots,\mathcal{F}_{K})$, and let the associated particle numbers be $(n_{1},n_{2},\dots,n_{K})$. The initial small-particle configuration is distributed as

[TABLE]

which is a normalised probability density. Then write

[TABLE]

where $q_{k}(\mathcal{F}_{k}\to\mathcal{F}_{k+1})$ is the probability (density) that the small particles evolve from $\mathcal{F}_{k}$ to $\mathcal{F}_{k+1}$, between steps $k$ and $k+1$. Rearranging factors in the exponential term, this becomes

[TABLE]

Detailed balance of the Markov chains implies that

[TABLE]

Using this in (65) together with (63) yields

[TABLE]

The integrals can now be performed one by one (starting with $\mathcal{F}_{1}$) and using that $q_{k}(\mathcal{F}_{k+1}\to\mathcal{F}_{k})$ is a normalised probability density for $\mathcal{F}_{k}$. The final integral over $\mathcal{F}_{K}$ is of the form (62) with $\mu=\mu_{\rm S}$ and one obtains (61).

Appendix B Simulation details

For simulation of the coarse model we use grand-canonical MC simulations as in Ref. Wilding, 1995. We restrict to particle insertion and removal (that is, no translational moves): this enables accurate results for fluid phases but helps to suppress crystallisation Largo and Wilding (2006). We perform long MC runs (for example up to $10^{8}$ MCS, for which the computational effort is a few hours) and we extract $\mathcal{N}_{\rm f}=128,000$ configurations for the computation of $\hat{\Delta}$.

The Jarzynski integration follows Sec. III.4 and Appendix A.2. The parameter $\mu_{0}$ is chosen such that

[TABLE]

with $\overline{n}_{0}$ of order unity being the average number of small particles that would be present, if there were no large particles. From grand-canonical simulations at this chemical potential (with fixed large particle configuration $\mathcal{C}$), we estimate $P_{\mathcal{C}}(n=0)$ as the probability to have no small particles at all, from which we compute $\Xi_{\rm S}(\mu_{0},\mathcal{C})=1/P_{\mathcal{C}}(n=0)$. The chemical potential is increased in steps with $\beta\mu_{k}=\beta\mu_{0}+\log(k+1)$, so that $(L/\sigma_{\rm S})^{3}\mathrm{e}^{\beta\mu_{k}}=(k+1)\overline{n}_{0}$. That is, the average number of particles in the system increases linearly with the number of steps. Between each step in $\mu$, we perform 1 MCS for the small particles. We take $\overline{n}_{0}=(0.25,0.48,0.45)$ for $q=(\frac{2}{5},\frac{1}{4},\frac{2}{13})$. The number of steps $K$ is fixed by the condition that $\mu_{K}=\mu_{\rm S}$, which is the chemical potential for the state point of interest.

The nature of the AO system means that accurate computation of the Jarzynski integral can be simplified. The system is divided into cells of size $\sigma_{\rm S}$. Before starting the computation, we identify cells that do not overlap at all with any large particle, and also cells where the large particles block the insertion of any small particle. We do not attempt to insert small particles in any of these cells. The remaining cells are called interaction cells and the number of them is $n_{\rm int}$. Hence the grand-canonical acceptance probabilities for particle insertion $P_{\rm in}$ and removal $P_{\rm del}$ are

[TABLE]

where $n^{\prime}$ is the number of small particles in the system. The time taken for the Jarzynski integration depends strongly on the number of small particles that have to be inserted, but a single calculation might take a few minutes, in a typical case.

We compute the Jarzysnki integral based on (24) with this MC method, which yields an unbiased estimate $\mathrm{e}^{\beta{\cal I}(\mathcal{C})}$ of $\exp[(\mathrm{e}^{\beta\mu_{\rm S}}-\mathrm{e}^{\beta\mu_{0}}){\cal V}^{\rm int}_{\rm a}(C)/\sigma_{\rm S}^{3}]$. Here, ${\cal V}^{\rm int}_{\rm a}$ is the volume accessible to the small particles, within the interaction cells. The total accessible volume is then ${\cal V}_{\rm a}=n_{\rm free}v_{\rm cell}+{\cal V}^{\rm int}_{\rm a}$, where $n_{\rm free}$ is the number of cells that do not overlap at all with any large particle. This number is known and so we can replace (25) with

[TABLE]

which still satisfies (19).

Appendix C Coarse model with three-body interactions

This Appendix gives extra information on three-body interactions in the system discussed in Sec. IV.3, which has $q=\frac{2}{5}$. As a starting point, Fig. 6 shows the performance of the TL method, using the coarse model of (14), which has only two-body interactions. In contrast to Fig. 4, it is clear that the TL method is not able to correct the coarse-graining errors in this coarse model. The failure of the method appears as a small number of very large weights $\hat{w}_{i}$ in Fig. 6c, up to $10^{10}$. The computational effort of Fig. 6 is comparable with that of Fig. 4, and $\mathcal{N}_{\rm f}=128,000$ in both cases. In principle the TL method can be used with any coarse model, but it is clear from (8) that $\mathcal{N}_{\rm f}$ should be significantly larger than the largest weights appearing in the calculation, otherwise large errors are expected. Hence one sees that the TL method is not feasible for the case shown in Fig. 6, since one would need $\mathcal{N}_{\rm f}\gtrsim 10^{10}$.

As discussed in Sec. IV.3, we analyse the system at $q=(2/5)$ using a CG model with (2+3)-body interactions. This interaction is represented as in (28), which is a scalar function $U_{3}$ of three scalar variables (which are interparticle distances). The function is symmetric in all its arguments and is non-zero only if all arguments are between $\sigma_{\rm B}$ and $\sigma_{\rm B}+\sigma_{\rm S}$. In practice we make the slightly stronger restriction that $U_{3}(x,y,z)=0$ except if $\sigma_{\rm B}^{2}<(x^{2}+y^{2}+z^{2})<(\sigma_{\rm B}+\sigma_{\rm S})^{2}$. We tabulate this function by estimating $U_{3}(x_{i},x_{j},x_{k})$ with $x_{j}=j\Delta\sigma$ and $\Delta\sigma=\sigma_{\rm S}/10$, and then we estimate the function at generic arguments using linear interpolation.

The tabulated values of the function are evaluated by Jarzynski integration as in (19), for systems with $N=3$ large particles. We take $L=3.6\sigma_{\rm B}$ and we generate $\mathcal{N}_{\rm f}=512,000$ configurations of the large particles by MC simulation of the coarse model (14). We perform Jarzynski integration on each of these configurations $\mathcal{C}$ which yields an estimate of $V_{3}(\mathcal{C})={6\eta_{\rm S}^{\rm r}}{\cal V}_{\rm a}(\mathcal{C})/(\pi\sigma_{\rm S}^{3}\beta)+U_{\rm c}(\mathcal{C})$ which is the three-body interaction in (17). For each $\mathcal{C}$ we compute the distances $R_{12},R_{13},R_{23}$ which are the arguments of $U_{3}$. Then we define $U_{3}(x_{i},x_{j},x_{k})$ as an average the estimated $V_{3}(\mathcal{C})$ over configurations $\mathcal{C}$ for which the 3-vector $(R_{12},R_{13},R_{23})$ is inside a cubic box of side $\Delta\sigma$, centred at $(x_{i},x_{j},x_{k})$.

An alternative approach would be to compute $U_{3}$ based on a deterministically chosen set of large-particle configurations. This might reduce the computational effort but it requires prior knowledge of the structure of the three-body potential. Instead, we computed the function at a set of representative points and we used these results to build the table. This “human-learning” method is sufficient for this simple three-body potential; machine-learning strategies would also be possible Behler and Parrinello (2007); Bartók et al. (2010); Cheng et al. (2019).

To illustrate the structure of the three-body potential, Fig. 7 shows the energy of the (2+3)-body coarse model for three particles in contact, as a function of the bond angle $\theta$. This calculation is performed by explicitly constructing large-particle configurations as in Fig. 7(c), and computing their weight factors using (19). Then $W_{\rm AO}(\mathcal{C})=\frac{6\eta_{\rm S}^{\rm r}}{\pi\sigma_{\rm S}^{3}}[{\cal V}_{\rm a}(\mathcal{C})-{\cal V}_{\rm a}(\mathcal{C}_{0})]$ where $\mathcal{C}_{0}$ is a configuration where all particles in the system are well-separated. This is the exact coarse energy of the (2+3)-body system, which is to be is compared with the energy of the two-body coarse model. Small $\theta$ (close to $\pi/3$) corresponds to the case where all three particles are close to each other. For $q=(1/2)$, the three-body interaction acts to increase the energy of such configurations, so the effective interaction between large particles is weaker than the prediction of the two-body coarse-model (14). For $q=(1/4)$ this effect is weaker, although the result of Fig. 3b shows that even these small changes in the interaction potential have significant effects close to criticality.

Appendix D The bias of the estimator $\hat{A}$

We show under the assumptions of Sec. V.1 that the expectation of $\hat{\mathcal{E}}$ appearing in (40) is of order $\mathcal{N}_{\rm f}^{-1}$. By definition, $\hat{\mathcal{E}}$ is an importance sampling estimator which converges to [math] almost surely for $\mathcal{N}_{\rm f}\to\infty$. In (38) we have seen that the normalisation of $\hat{w}$ slightly perturbs the weights appearing in $\hat{\mathcal{E}}$ , i.e. in general $\epsilon\neq 0$. This results in a bias which we can estimate by following Ref. Agapiou et al., 2017. For completeness, we reproduce here the steps used in the proof of their Theorem 2.1, adapted to our notation.

We introduce the short-hand notation

[TABLE]

Also note from (19) that

[TABLE]

Using the definition of $\hat{\cal E}$ with (36,38,73) yields

[TABLE]

We decompose the expectation value into contributions from $\epsilon\leq-\frac{1}{2}$ and $\epsilon>-\frac{1}{2}$, so that

[TABLE]

where the indicator function $1_{\{X\}}$ is equal to unity if $X$ is true, and zero otherwise. For the first term on the right hand side of (75), we use that $|\epsilon/(1+\epsilon)|\leq 2|\epsilon|$ for $\epsilon>-\frac{1}{2}$ to obtain

[TABLE]

The right hand side is less than or equal to

[TABLE]

which follows from the Cauchy-Schwarz inequality. Using this fact together with (37) and that the $\mathcal{C}_{i}$ are all independent, we find

[TABLE]

Then (34) yields

[TABLE]

For the second term on the right hand side of (75), one uses (20,34,41) to show that $|\hat{\mathcal{E}}|\leq a$. Hence

[TABLE]

This probability can be bounded by Chebyshev’s inequality which implies that for any mean-zero random variable $X$ and any number $x>0$ then $\mathrm{Prob}(X\leq-x)\leq\mathrm{Var}(X)/x^{2}$. Hence

[TABLE]

Combining (75,79,81) we have

[TABLE]

which is an explicit bound that holds for all $\mathcal{N}_{\rm f}$. The right hand side is finite by (33) so $\langle\hat{\mathcal{E}}\rangle$ is $O(\mathcal{N}_{\rm f}^{-1})$ as advertised above. With (40) this establishes (42).

We note in passing that if the weights have zero variance then $\hat{\cal E}=0$ strictly and there is no bias, although the bound (82) does not capture this fact. In fact this bound may be improved by refining the estimate of $\langle\epsilon^{2}\rangle$ used in (78).