Thermodynamics of mixtures of patchy and spherical colloids of different sizes: a multi-body association theory with complete reference fluid information

TL;DR

This paper develops a multi-body association theory for predicting the structure and thermodynamics of mixtures of patchy and spherical colloids of different sizes, incorporating complete reference fluid information for improved accuracy.

Contribution

It introduces a comprehensive multi-body association framework that accounts for full reference fluid correlations, enhancing predictions of colloid mixture behavior over previous models.

Findings

Complete reference information improves thermodynamic predictions.

The theory accurately predicts bonding states across various conditions.

Enhanced sampling and statistical models estimate occupancy probabilities effectively.

Abstract

We present a theory to predict the structure and thermodynamics of mixtures of colloids of different diameters, building on our earlier work [J. Chem. Phys. 145, 074904 (2016)] that considered mixtures with all particles constrained to have the same size. The patchy, solvent particles have short-range directional interactions, while the solute particles have short-range isotropic interactions. The hard-sphere mixture without any association site forms the reference fluid. An important ingredient within the multi-body association theory is the description of clustering of the reference solvent around the reference solute. Here we account for the physical, multi-body clusters of the reference solvent around the reference solute in terms of occupancy statistics in a defined observation volume. These occupancy probabilities are obtained from enhanced sampling simulations, but we also present statistical mechanical models to estimate these probabilities with limited simulation data. Relative to an approach that describes only up to three-body correlations in the reference, incorporating the complete reference information better predicts the bonding state and thermodynamics of the physical solute for a wide range of system conditions. Importantly, analysis of the residual chemical potential of the infinitely dilute solute from molecular simulation and theory shows that whereas the chemical potential is somewhat insensitive to the description of the structure of the reference fluid the energetic and entropic contributions are not, with the results from the complete reference approach being in better agreement with particle simulations.

Tables

Table 1. Table 1: Parameters for Eq.  24 ( ζ 1 subscript 𝜁 1 \zeta_{1} , ζ 2 subscript 𝜁 2 \zeta_{2} ) and Eq.  26 ( ζ 𝜁 \zeta ) for different densities ρ ​ σ 3 𝜌 superscript 𝜎 3 \rho\sigma^{3} and packing fractions ( η 𝜂 \eta )

$\rho {\sigma }^3$	$\eta$	${\zeta }_1$	${\zeta }_2$	$\zeta$
0.2	0.105	0.0175	0.773	0.75
0.3	0.157	0.015	1.267	1.251
0.4	0.209	0.0256	1.979	1.947
0.5	0.262	0.023	2.88	2.876
0.6	0.314	0.0361	4.179	4.172
0.7	0.367	0.0457	5.947	5.985
0.8	0.419	0.0609	8.432	8.562
0.9	0.471	0.0829	12.088	12.414

Table 2. Table 2: Comparison of average bonding numbers ( n a ​ v ​ g subscript 𝑛 𝑎 𝑣 𝑔 n_{avg} ) of solute for extreme size ratios.

${\sigma }_s/{\sigma }_p$	MC	$p_n$-Simulation	$p_n$-Model1
0.5	3.61	3.55	4.50
2	15.84	15.04	17.66

Table 3. Figure 11: Left: Maximum number of coordination states possible and their respective maximum permissible patch angle for different size ratios σ s / σ p subscript 𝜎 𝑠 subscript 𝜎 𝑝 \sigma_{s}/\sigma_{p} . Right: Diagram to compute the angle between two solvent particles to get the maximum number of coordination states

${\sigma }_s/{\sigma }_p$	$r_c$	$n^{max}$	${\theta }_{c,max}$
0.5	0.80	8	17.64
0.8	0.99	12	23.83
0.9	1.05	12	25.51
1.0	1.10	14	27.04
1.1	1.16	16	28.44
1.2	1.21	18	29.73
2.0	1.55	31	40.18

Table 4. Table 3: Isolated cluster probabilities P ( n ) superscript 𝑃 𝑛 P^{(n)} for different size ratios.

	${\sigma }_s/{\sigma }_p$
$n$	0.8	0.85	0.9	0.95	1	1.05	1.1	1.15	1.2
$1$	$1$	$1$	$1$	$1$	$1$	$1$	$1$	$1$	$1$
$2$	$0.72087$	$0.73574$	$0.74949$	$0.76219$	$0.77396$	$0.78485$	$0.79497$	$0.80444$	$0.81326$
$3$	$0.34732$	$0.37301$	$0.39762$	$0.42114$	$0.4435$	$0.46487$	$0.48523$	$0.50462$	$0.52305$
$4$	$0.10073$	$0.11924$	$0.13835$	$0.15792$	$0.1778$	$0.19765$	$0.21762$	$0.23745$	$0.25702$
$5$	$0.01518$	$0.02126$	$0.02854$	$0.03689$	$0.04632$	$0.05668$	$0.06793$	$0.07992$	$0.0926$
$6$	$9.61\cdot {10}^{-4}$	$1.78\cdot {10}^{-3}$	$3.03\cdot {10}^{-3}$	$4.77\cdot {10}^{-3}$	$7.1\cdot {10}^{-3}$	$0.01008$	$0.01374$	$0.01812$	$0.02322$
$7$	$1.73\cdot {10}^{-5}$	$5.28\cdot {10}^{-5}$	$1.33\cdot {10}^{-4}$	$2.89\cdot {10}^{-4}$	$5.61\cdot {10}^{-4}$	$9.97\cdot {10}^{-4}$	$1.64\cdot {10}^{-3}$	$2.56\cdot {10}^{-3}$	$3.79\cdot {10}^{-3}$
$8$	$6.04\cdot {10}^{-8}$	$4.01\cdot {10}^{-7}$	$1.86\cdot {10}^{-6}$	$6.59\cdot {10}^{-6}$	$1.91\cdot {10}^{-5}$	$4.72\cdot {10}^{-5}$	$1.03\cdot {10}^{-4}$	$2.03\cdot {10}^{-4}$	$3.7\cdot {10}^{-4}$
$9$	$1.21\cdot {10}^{-10}$	$4.25\cdot {10}^{-10}$	$5.26\cdot {10}^{-9}$	$3.98\cdot {10}^{-8}$	$2.12\cdot {10}^{-7}$	$8.64\cdot {10}^{-7}$	$2.86\cdot {10}^{-6}$	$7.95\cdot {10}^{-6}$	$1.93\cdot {10}^{-5}$
$10$	$1.17\cdot {10}^{-14}$	$3.77\cdot {10}^{-13}$	$3.56\cdot {10}^{-12}$	$8.59\cdot {10}^{-11}$	$1.48\cdot {10}^{-9}$	$4.37\cdot {10}^{-9}$	$2.74\cdot {10}^{-8}$	$1.27\cdot {10}^{-7}$	$4.67\cdot {10}^{-7}$
$11$		$1.73\cdot {10}^{-17}$	$2.57\cdot {10}^{-15}$	$5.5\cdot {10}^{-14}$	$4.98\cdot {10}^{-12}$	$3.87\cdot {10}^{-12}$	$6.62\cdot {10}^{-11}$	$6.46\cdot {10}^{-10}$	$4.3\cdot {10}^{-9}$
$12$		$6.33\cdot {10}^{-24}$	$2.6\cdot {10}^{-20}$	$5.57\cdot {10}^{-18}$	$1.11\cdot {10}^{-15}$	$1.42\cdot {10}^{-15}$	$4.78\cdot {10}^{-14}$	$2.29\cdot {10}^{-12}$	$1.27\cdot {10}^{-11}$
$13$					$5.55\cdot {10}^{-24}$	$4.44\cdot {10}^{-19}$	$4.27\cdot {10}^{-17}$	$1.07\cdot {10}^{-15}$	$4.77\cdot {10}^{-15}$

Equations

$$u(r)=u_{HS}(r)+u_{AS}(r)$$

$$u_{AB}^{(p,p)}(r)=\{\begin{array}{l}-{ϵ}_{AB}^{(p,p)},r<r_c\text{and}{\theta }_A\le {\theta }_c^{(A)}\text{and}{\theta }_B\le {\theta }_c^{(B)}\\ 0\text{ otherwise}\end{array}$$

$$u_A^{(s,p)}(r)=\{\begin{array}{l}-{ϵ}_A^{(s,p)},r<r_c\text{and}{\theta }_A\le {\theta }_c^{(A)}\\ 0\text{ otherwise}\end{array}$$

$$F^{(n)}$$

$$\frac{A^{AS}}{V{k}_{\mathrm{B}}T}=\sum \left({\rho }_k\ln \frac{{\rho }_k^{(0)}}{{\rho }_k}+Q^{\left(k\right)}+{\rho }_k\right)-\frac{\Delta c^{(0)}}{V}$$

$$\Delta c^{\left(0\right)}=\Delta c_{pp}^{\left(0\right)}+\Delta c_{sp}^{\left(0\right)}$$

$$\frac{\Delta c_{sp}^{(0)}}{V}=\sum _{n=1}^{n^{\max }}\frac{\Delta c_n^{(0)}}{V}$$

$$\Delta c_n^{(0)}$$

$$\frac{\Delta c_n^{(0)}}{V}$$

$${\int }_{v}d{\vec{r}}_1\cdots {\int }_{v}d{\vec{r}}_n{g}_{HS}({\vec{r}}_1\cdots {\vec{r}}_n|0)\approx y_{HS}^n(\sigma ){\delta }^{(n)}{\Xi }^{(n)},$$

$$\frac{\Delta c_n^{(0)}}{V}={\rho }_s^{(0)}{(x^{(p)}{X}_A^{(p)}{f}_A^{(s,p)}\sqrt{{\kappa }_{AA}})}^n{F}^{(n)}.$$

$$X_n^{\left(s\right)}=\frac{{(x^{(p)}{X}_A^{(p)}{f}_A^{(s,p)}\sqrt{{\kappa }_{AA}})}^n{F}^{(n)}}{1+\sum _{n=1}^{n^{\max }}{(x^{(p)}{X}_A^{(p)}{f}_A^{(s,p)}\sqrt{{\kappa }_{AA}})}^n{F}^{(n)}},$$

$$X_0^{\left(s\right)}=\frac{1}{1+\sum _{n=1}^{n^{\max }}{(x^{(p)}{X}_A^{(p)}{f}_A^{(s,p)}\sqrt{{\kappa }_{AA}})}^n{F}^{(n)}}.$$

$$n_{avg}=\sum _{n}n\cdot X_n^{(s)},$$

$$X_A^{\left(p\right)}=\frac{1}{1+\xi {\kappa }_{AB}{f}_{AB}^{\left(p,p\right)}{\rho }_p{X}_B^{(p)}+\frac{{\rho }_s}{{\rho }_p}\frac{n_{avg}}{X_A^{(p)}}},$$

$$X_B^{\left(p\right)}=\frac{1}{1+\xi {\kappa }_{AB}{f}_{AB}^{\left(p,p\right)}{\rho }_p{X}_A^{(p)}}.$$

$$\xi$$

$${\kappa }_{AB}$$

$$f_{AB}^{(p,p)}$$

$$S{P}_{n=0}+P_n\rightleftharpoons S{P}_n,$$

$$K_n=\frac{{\rho }_{S{P}_n}}{{\rho }_{S{P}_{n=0}}{\rho }_p^n},$$

$$\frac{p_n}{p_0}=K_n{\rho }_p^n,$$

$$K_n=\frac{(e^{\beta {\mu }_p^{\mathrm{ex}}}{)}^n}{n!}\underset{v}{\int }d{\vec{r}}_1\dots \underset{v}{\int }d{\vec{r}}_n{e}^{-\beta U_{S{P}_n}(R^n)}{e}^{-\beta \varphi (R^n;\beta )},$$

$$K_n=\frac{(e^{\beta {\mu }_p^{\mathrm{ex}}}{)}^n}{n!}⟨e^{-\beta \varphi (R^n;\beta )}{⟩}_0\underset{v}{\int }d{\vec{r}}_1\dots \underset{v}{\int }d{\vec{r}}_n{e}^{-\beta U_{S{P}_n}(R^n)}.$$

$$n!K_n^{(0)}=\underset{v}{\int }d{\vec{r}}_1\dots \underset{v}{\int }d{\vec{r}}_n{e}^{-\beta U_{S{P}_n}(R^n)},$$

$$\beta {\Omega }_n=-{\zeta }_1\cdot n^2+{\zeta }_2\cdot n.$$

$${p_{n}}=\frac{\exp(\beta\cdot n\cdot\mu^{ex}_{p})\rho_{p}^{n}[{\exp{(\zeta_{1}\cdot n^{2}-\zeta_{2}\cdot n)}}]{K_{n}}^{(0)}}{{1+\sum\limits_{j\geq 1}{\exp(\beta\cdot j\cdot\mu^{ex}_{p})\rho_{p}^{j}[{\exp{(\zeta_{1}\cdot j^{2}-\zeta_{2}\cdot j)}}]{K_{j}}^{(0)}}}}\$$

$$\beta {\Omega }_n=(-0.0109\cdot n^2+1.0109\cdot n)\cdot \zeta$$

$${p_{n}}=\frac{\exp(\beta\cdot n\cdot\mu^{ex}_{p})\rho_{p}^{n}[{\exp{(-(-0.0109\cdot n^{2}+1.0109\cdot n)\cdot\zeta)}}]{K_{n}}^{(0)}}{{1+\sum\limits_{j\geq 1}{\exp(\beta\cdot j\cdot\mu^{ex}_{p})\rho_{p}^{j}[{\exp{(-(-0.0109\cdot j^{2}+1.0109\cdot j)\cdot\zeta)}}]{K_{j}}^{(0)}}}}\$$

$$\beta {\mu }_s^{Asso}=\beta E_s^{Asso}-T\cdot \beta S_{s,V}^{Asso}$$

$$\frac{p_n}{p_0}=\frac{N!}{(N-n)!n!}\frac{\underset{v}{\int }d{\vec{r}}_1\dots \underset{v}{\int }d{\vec{r}}_n{e}^{-\beta U_{S{P}_n}}\cdot \underset{V-v}{\int }d{\vec{r}}_{n+1}\dots \underset{V-v}{\int }d{\vec{r}}_{N-n}{e}^{-\beta U_{N-n|S{P}_n}}{e}^{-\beta U_{N-n}}}{\underset{V-v}{\int }d{\vec{r}}_1\dots \underset{V-v}{\int }d{\vec{r}}_N{e}^{-\beta U_N}}$$

Full text

I Introduction

In thermodynamic perturbation theory of association involving short range interactions between molecules, the properties of the reference fluid plays a central role. In the typical situation when the reference is a hard-sphere fluid, perturbation theories usually use information about two body, and at times three body, correlations in the reference fluid to describe the physical (associating) system. For example, Wertheim’s theory Wertheim (1984a, b) and its extensions based on the statistical associating fluid theory (SAFT) Chapman et al. (1990) use pair correlation information at contact to estimate extent of association between pairs of molecules. In SAFT, for the hard-sphere reference either the Carnahan-Starling Carnahan and Starling (1969) equation for a single component fluid or the Boublik-Mansoori-Carnahan-Starling-Leland Boublík (1970); Mansoori et al. (1971) equation for a mixture are used to describe the pair-correlation information at contact. The structure of Wertheim’s theory or SAFT is such that one can obtain accurate extent of association and thermodynamics even for systems with strong inter-particle interactions provided the representation of the reference is adequate. However, as the complexity of the interaction increases in the physical system, such as may result from multiple bonding and size asymmetries, information about two or three body correlations in the reference no longer suffices.

In our previous work Bansal et al. (2016), we studied the multi-body correlation functions of a symmetric hard sphere reference fluid in terms of the probabilities of observing $n$ molecules in the bonding region. These occupancy probabilities were obtained from enhanced sampling Monte Carlo simulations for the hard sphere fluid. We developed a procedure to use this information within the Marshall-Chapman formalism Marshall and Chapman (2013a, b) to describe multiple association of solvent molecules to a solute molecule. This complete reference approach proved successful in predicting the bonding state and thermodynamics of a colloidal solute in a patchy solvent for a wide range of system conditions Bansal et al. (2016).

Here we study mixtures where the solute diameter is as small as half to as large as twice the diameter of the solvent. The solvent particles are spheres with directional interaction sites and the solute particles are spheres with isotropic interactions, and the solute is capable of bonding with multiple solvent particles. The structure and thermodynamics of mixture of hard spheres with different diameters has been studied in detail before Torquato (1986); Reiss et al. (1959); Mayer (1947); Torquato and Stell (1982, 1983, 1985), but a compact form for the correlations beyond the contact value is still unavailable. Further, for systems with large size asymmetries even the pair-correlation information obtained using the Boublik-Mansoori-Carnahan-Starling-Leland equation is inadequate Feng and Chapman (2011). Our approach of including multi-body correlations rests on using the occupancy statistics Reiss and Merry (1981); Pratt et al. (2001); Pratt and Ashbaugh (2003); Bansal et al. (2016) of the hard-sphere solvent around the hard-sphere solute. We find that representing multi-body correlation functions in terms of occupancy statistics in physically reasonable observation volumes accurately captures the multi-bonding effects in asymmetric mixtures. These occupancy statistics are obtained from particle simulations. Importantly, we also present a physically transparent, statistical mechanical model to describe the occupancy probabilities in symmetric and asymmetric hard sphere fluids for different packing fractions. This model corrects multi-body effects obtained for isolated clusters by incorporating the role of the cluster-bulk interface and the bulk medium effects. We also investigate the energy-entropy decomposition of the chemical potential of the solute in a model system with only solute-solvent interactions to better appreciate the role of the reference fluid. Throughout, theoretical results are validated versus molecular simulations

The rest of the paper is organized in the following way. In Section II.1 we discuss the association potential of the system and describe how packing effects are important for the given potential. The Marshall-Chapman Marshall and Chapman (2013b) theory is briefly introduced to show the multi-density representation of the free energy, and based on our previous work Bansal et al. (2016), an improved representation of multi-body correlations(complete reference) is presented. In Section II.2, we examine hard sphere packing around a reference particle and develop models based on statistical mechanics Reiss and Merry (1981) and hard sphere simulation data for different densities. We apply our complete reference approach for different asymmetric mixtures of solute and solvent and present results in Section IV.2. In Section IV.1 we present results for the hard sphere reference system (symmetric and asymmetric mixtures) based on the correlation developed in section II.2. We also provide simulation results for isolated cluster probabilities in asymmetric hard sphere mixtures in the appendix (Section VII).

II Theory

II.1 Asymmetric mixtures with different association geometries

The focus of our study is asymmetric mixtures containing molecules with short range attractive interactions. The short range association potential is the same as that in previous work Bansal et al. (2016): the solute molecule can associate with multiple solvent molecules isotropically and the patchy solvent has directional interactions. The total potential is a sum of hard sphere and association contributions

[TABLE]

The association potential for patchy-patchy $(p,p)$ and spherical-patchy $(s,p)$ particles is:

[TABLE]

[TABLE]

where the subscripts $A$ and $B$ represent the type of site and $\epsilon$ is the association energy; $r$ is the distance between the particles; and $\theta_{A}$ is the angle between the vector connecting the centers of two molecules and the vector connecting association site $A$ to the center of that molecule (Fig. 1). The critical distance beyond which particles do not interact is $r_{c}$ and $\theta_{c}$ is the solid angle beyond which sites cannot bond. Fig. 1 shows examples of solute-solvent and solvent-solvent short range interaction geometries for different sizes of solute particles.

Since the solute can associate with multiple solvent molecules (Eq. 3), it is important to study the multi-body correlations that determine the packing of solvent particles around the solute in the reference fluid Bansal et al. (2016). The difficulty in determining these interactions arises due to the limited knowledge in describing multi-body correlation functions for $n\geq 3$. But the volume integral of the multi-body correlation has a clear physical meaning in terms of average number of $n$-solvent clusters ($F^{(n)}$, Fig. 2). In particular, for the distinguished solute,

[TABLE]

where $\rho_{p}$ is the density of solvent particles, $p_{n}$ is the probability of observing exactly $n$ solvent particles in the observation volume of the solute ($v$) defined by the spherical region of radius $r_{c}$, $C^{m}_{n}$ ($=m!/(m-n)!\cdot n!$), and $g_{HS}({\vec{r}}_{1}\cdots{\vec{r}}_{n}|0)$ is the distribution function of the $n$-solvent particles around the solute at the center of the observation volume, indicated by $(\ldots|0)$. $n^{max}$ is the maximum number of solvent molecules that can occupy the observation volume around the reference solute.

In Wertheim’s multi-density formalism Wertheim (1986a, b), the free energy due to association ($A^{AS}$) is expressed as

[TABLE]

where $k_{\rm B}$ is the Boltzmann constant, $T$ is the temperature, the summation is over the species ($k={s,p}$), $\rho$ is the number density, $\rho^{(0)}$ is the monomer density, $Q^{(k)}$ is obtained from the Marshall-Chapman development Marshall and Chapman (2013b) and $\Delta{c^{0}}$ is the contribution to the graph sum due to association between the solvent-solvent $(p,p)$ and solute-solvent $(s,p)$ molecules, i.e.

[TABLE]

Marshall and Chapman Marshall and Chapman (2013a, b) extended Wertheim’s theory beyond the single bonding condition to incorporate multi-body effects in a solution consisting of an isotropic solute and solvent with directional interactions. The contribution to free energy due to association between solute and solvent molecules was obtained as

[TABLE]

where the sum is over different coordination states of the solute and $\Delta c_{n}^{(0)}$ is given by:

[TABLE]

In Eq. 8, $\rho_{p}=\rho\cdot x^{(p)}$ is the density of solvent molecules obtained from the mole fraction of solvent($x^{(p)}$) and the total density($\rho$), $X_{A}^{(p)}$ is the fraction of solvent molecules not bonded at site A, $\Omega=4\pi$ is the total number of orientations, $f_{A}^{({s,p})}(1,k)=(\exp(\varepsilon_{A}^{(s,p)}/k_{B}T)-1)$ is the Mayer function for association between $p$ and $s$ molecules corresponding to potential in Eq. 3 and the integral is over all the orientations and positions of the $n+1$ particles. By taking the average association strength and acceptable orientations out of the integral and fixing the solute at the origin, the above integral can be rewritten as

[TABLE]

Marshall and ChapmanMarshall and Chapman (2013a, b) approximated the integral in Eq. 9 as

[TABLE]

where $\Xi^{(n)}$ is the partition function for an isolated cluster of $n$ solvent hard-spheres around a solute hard-sphere, $y_{HS}(\sigma)$ is (pair) cavity correlation function at contact, and $\delta^{\left(n\right)}$ corrects the superposition of cavity correlation functions for three body interactions. We will hereafter refer to Eq. 10 as the Marshall-Chapman approximation (MCA).

As shown earlier Bansal et al. (2016), MCA fails for high densities and high association energies, conditions where multi-body interactions are important. But recognizing that the integral in Eq. 9 is related to $F^{(n)}$ (Eq. 4) we have Bansal et al. (2016)

[TABLE]

It can be observed that all the multi-body correlation information is subsumed in $F^{(n)}$ which is obtained from the occupancy distribution $\{p_{n}\}$. We follow our earlier work Bansal et al. (2016) to estimate this distribution. Importantly, since $\{p_{n}\}$ forms the basis of our complete reference approach, we also develop an analytical model to describe these distribution functions.

Finally, with the above information, and based on the Marshall-Chapman theory Marshall and Chapman (2013b), the fraction of solute associated with $n$ solvent molecules is

[TABLE]

and the fraction of solute not bonded to any solvent molecule is

[TABLE]

Using these distributions for associating mixture, the average number of solvent associated with the solute is given by:

[TABLE]

The fraction of solvent not bonded at site $A$ and site $B$ can be obtained by simultaneous solution of the following equations:

[TABLE]

[TABLE]

where

[TABLE]

II.2 Occupancy distribution $\{p_{n}\}$ for the hard-sphere fluid

Consider a hard sphere fluid with one solute and $N$ solvent particles in a volume $V$ and temperature $T$. We are interested in the occupancy statistics $\{p_{n}\}$ of the solvent in the coordination volume around the solute. To this end consider the reaction

[TABLE]

with the equilibrium constant

[TABLE]

where $\rho_{SP_{n}}$ is the density of species $SP_{n}$ and $\rho_{p}$ is the density of the solvent. Clearly, we have Pratt et al. (2001); Beck et al. (2006); Pratt and Asthagiri (2007)

[TABLE]

where $p_{0}$ is the probability that the coordination volume is empty of solvent particles. Following earlier work in studying clusters with quasichemical theory Pratt et al. (2001); Pratt and Ashbaugh (2003); Merchant and Asthagiri (2009); Merchant et al. (2011a), we can show that

[TABLE]

where $U_{SP_{n}}(R^{n})$ is the potential energy of the $n$-solvent cluster (with the solute $S$ fixed at the center of the cluster), $\beta=(k_{\rm B}T)^{-1}$, $\phi(R^{n};\beta)$ is the free energy of interaction of the cluster with the rest of the bulk medium for a given configuration $R^{n}$ of the cluster, and $\beta\mu^{\rm ex}_{p}$ is the excess chemical potential of the solvent particle. (For completeness, in appendix A we derive the above expression for $K_{n}$.) $\phi(R^{n};\beta)$ can also be thought of as a field imposed by the bulk solvent medium Pratt and Ashbaugh (2003); Merchant et al. (2011a) on the solute-solvent cluster in the observation volume. Earlier Pratt and Ashbaugh Pratt and Ashbaugh (2003) modeled this field using a self-consistent approach. Here we take a different approach.

First note that without the field term, the cluster integral presents a simpler $n$-body problem (where $n$ is small, typically less than 20 for systems of interest here). The field is thus an interfacial term that couples the local cluster with the bulk medium. To make this explicit, we can rewrite Eq. 20 as

[TABLE]

Here $\langle\ldots\rangle_{0}$ indicates averaging with over the normalized probability density for cluster conformations $R^{n}$ in the absence of interactions with the rest of the medium, i.e. over the density $e^{-\beta U_{SP_{n}}(R^{n})}/(n!K_{n}^{(0)})$, where

[TABLE]

and the interfacial contribution is $\beta\Omega_{n}=-\ln\langle e^{-\beta\phi(R^{n};\beta)}\rangle_{0}$.

From analysis of simulation data for different densities, we find that $\Omega_{n}$ can be described by a two parameter equation as

[TABLE]

This model of the interfacial term was anticipated in our previous work Bansal et al. (2016). Thus we finally obtain

[TABLE]

The above equation can also be derived using a two moment maximum entropy approach, with the mean and variance of the occupancy as constraints and $K_{n}^{(0)}$ as the default (see appendix A).

Drawing upon the work of Reiss and Merry Reiss and Merry (1981), we can model the interfacial term in terms of surface sites (of the cluster) that are available to interact with the bulk fluid. On the basis of such a mean field approach and guided by Monte Carlo (MC) simulation data for different densities of hard sphere systems, we find that

[TABLE]

[TABLE]

Eq. 24 and Eq. 26 are the 2-parameter and 1-parameter models, respectively, for $p_{n}$ on the basis of which we obtain $F^{(n)}$ (Eq. 4) to describe multi-body correlations in the reference fluid. The parameter values for different densities are given in Table 1. These parameters were obtained based on hard spheres mixtures with all particles of the same size. Since the information about size-asymmetry is already contained in the isolated cluster partition function, we will use these parameters to study asymmetric mixtures and will discuss limitations for cases with extreme size ratio.

III METHODS

To compare the theory results, we perform Monte-Carlo (MC) simulations for a range of systems. This section presents the details of the MC simulations for different associating and hard sphere systems.The Marshall-Chapman approximation (MCA) and the models developed for hard sphere distribution functions require isolated cluster probabilities; these were also computed for different size ratios.

III.1 Monte Carlo Simulations

MC simulations were carried out for reference hard sphere systems and associating systems to compare the results of Marshall-Chapman theory using MCA and the complete reference approach. The associating mixture contains the patchy solvent particles and the isotropically interacting solute defined by the potentials given by Eq. 2 and Eq. 3, respectively. The solute diameter is $\sigma_{s}$ and solvent diameter is $\sigma_{p}$.

The observation volume is defined by a critical radius $r_{c}=1.1\bar{\sigma}$, where $\bar{\sigma}=(\sigma_{s}+\sigma_{p})/2$ is the closest distance of approach. For cases where $\sigma_{s}/\sigma_{p}\geq 1.5$, a cutoff of 1.1$\bar{\sigma}$ can include some of the second-shell solvent. To avoid this and focus attention to the first observation shell, we set $r_{c}=\bar{\sigma}+0.1\sigma_{p}$ for these cases. In the associating system (Fig. 1 and Eqs. 2 and 3), the critical angles for interaction are $\theta_{c}^{(A)}=\theta_{c}^{(B)}=27\degree$.

For hard sphere mixtures 255 solvent particles and 1 solute particle were studied in a given simulation cell. Ensemble reweighting technique was used to reveal low probability states Merchant et al. (2011b). The system was equilibrated for 1 million steps with translational factors chosen to yield an acceptance rate of 0.3, and data was collected every 100 sweeps. Analysis was carried out for different densities and size ratios $\sigma_{s}/\sigma_{p}$.

For associating mixtures bonding distributions and average bonding numbers were studied for mixtures with different sizes and different association strengths for solute-solvent and solvent-solvent interactions. The excess chemical potential of the coupling of the colloid with solvent was also calculated using thermodynamic integration of average binding energy of solute with solvent as a function of solute-solvent interactions, using the three-point Gauss Legendre quadrature technique Hummer and Szabo (1996). For a symmetric mixture with no solvent-solvent interactions, energetic and entropic contributions for solute chemical potential were also studied at constant volume and temperature.

Concentration effects were also computed considering a total of 864 particles, with varied number of solute particles. Due to the difference in size of the solute and solvent, the computations were performed keeping the packing fraction constant. Hence, the density of the system changed with respect to the concentration. Also the maximum angle for which the patch can form single bond is computed using the law of cosines, and hence the critical angle needs to be altered when the solute size is smaller than the solvent. For a size ratio of $\sigma_{s}/\sigma_{p}=0.8$, critical angle $\theta_{c}^{(A)}=\theta_{c}^{(B)}=20\degree$ was used to ensure single bonding condition for the $A$ patch on solvent molecules which can associate with solute molecules.

III.2 Isolated cluster probabilities

For asymmetric mixtures, we also study isolated cluster probabilities for different size ratios of solute and solvent molecules. The observation volume around the isolated solute is the same as defined in the previous section. For different size ratios we use the spherical code (appendix B) to estimate the maximum number of solvent molecules which can be inserted in the observation volume of the isolated solute. Successive insertion probabilities are calculated as in previous works Marshall and Chapman (2013a); Bansal et al. (2016), where $10^{8}$ to $10^{9}$ insertions were carried out in the observation volume around the solute and the cases with no overlap with remaining $n-1$ particles studied. The data for isolated cluster probabilities for different sizes is given in the appendix B.

IV Results and discussions

IV.1 Hard sphere $\{p_{n}\}$ distribution

Recall that Eq. 26 and Eq. 24 are the 1-parameter and 2-parameter models for the occupancy distribution $\{p_{n}\}$. To obtain the parameters for hard spheres all of the same size, for both the models we use the average occupancy ($n^{HS}_{avg}=\sum\limits_{n}n\cdot{p_{n}}$) as a fitting constraint. For the 2-parameter model we additionally use the exclusion probability ($p_{0}$) — the probability when no hard sphere solvent particle is present in the observation volume — as a constraint. For the 1-parameter model we study the surface interactions based on the mean field approach developed by Reiss and Merry Reiss and Merry (1981). By analyzing the distribution functions $\{p_{n}\}$ for different densities, we obtain geometric effects (density independent) that describe the mutual interference of different surface sites (Eq. 25). The density (or packing fraction) dependent parameters for these two models are given in Table  1.

Fig. 3 presents the results corresponding to the average occupancy ($n^{HS}_{avg}$) and exclusion probability ($\ln p_{0}$) based on the models. We compare these results with the MC simulation values and also include results from literature Chang and Sandler (1994); Torquato and Stell (1985).

Fig. 3 makes it clear that 2-parameter model can simultaneously capture both $n_{avg}$ and $\ln p_{0}$ in excellent agreement with simulation. Importantly, even the 1-parameter model is able to capture most of the details, affirming the physical ideas underlying the models (Eqs. 26 and  24).

Next, using the parameters obtained above for a fluid where both solute and solvent hard-spheres are the same size (a symmetric mixture), we describe the occupancy in a fluid where the solute and solvent are of different sizes (an asymmetric mixture). Our ansatz is that the information about size asymmetry is adequately captured by the isolated cluster partition function $K_{n}^{(0)}$ (Eq. 22). For an infinitely dilute system comprising one solute in a solvent bath, Fig. 4 shows the predictions of $n_{avg}^{HS}$ and $\ln p_{0}$ based on the 1-parameter model (Eq. 26) for different size ratios and different reduced densities. The level of agreement is encouraging, but perhaps not surprising since the size ratios are not much different from a symmetric case. Thus the geometric effects describing the mutual interference of different surface sites in the packing around the solute will be similar to what is observed in the symmetric mixture. For extreme size ratios this should break down, as we will discuss in the context of average bonding numbers in associating mixtures.

For asymmetric mixtures, the packing fraction is a better measure of packing in the fluid. Fig. 5

presents the variation of average occupancy for different packing fractions for three different size ratios for an infinitely dilute solution. The results show that the 1-parameter model is able to predict the average occupancy quite well. We note that as the concentration of solute is changed in asymmetric mixtures, the packing fraction changes (for a given density) and parameters corresponding to the resulting packing fraction should be used from Table 1.

IV.2 Associating mixture

We next consider associating fluids and investigate both size and concentration effects.

IV.2.1 Infinite Dilution

We first study an infinitely dilute solution and vary the size of solute with respect to a fixed size of solvent particles. In our complete reference theory, the reference fluid $\{p_{n}\}$ distribution is either computed directly from simulations (‘$p_{n}$-Simulation’ in Fig. 6) or from the 1-parameter model discussed in Sec. A above (‘$p_{n}$-Model1’ in Fig. 6). Using $\{p_{n}\}$ we compute $F^{(n)}$ (Eq. 4), and on that basis, the average number of bonds in the associating mixture using Eq. 11, Eq. 12 and Eq. 14. Fig. 6 shows the variation of average bonding numbers with size ratio of solute and solvent molecules for a density of 0.8 and association strength of 7 $k_{\rm B}T$.

As the size of the solute increases with respect to the size of the solvent, more solvent molecules can associate with the solute. With accurate information about the reference fluid $\{p_{n}\}$ (and hence $F^{(n)}$) from direct simulations, the complete reference theory is able to capture this increase in bonding numbers quite accurately. Interestingly, even $\{p_{n}\}$ obtained using the 1-parameter model (Eq. 26) suffices. But note that the Marshall-Chapman approximation with only up to 3-body effects incorporated in the theory underestimates the average bonding numbers. These results emphasize the importance of multi-body interactions in describing the association correctly.

As noted in Sec.  II.2, the amount of surface exposure (or surface sites) is an important factor in determining the packing effects in the models. For the size ratios considered in Fig. 6, the maximum number of surface sites and hence the geometric effects in surface interactions are expected to be similar to the symmetric case, and not surprisingly, the agreement of bonding numbers between simulation and with the 1-parameter model for $p_{n}$ is very good. For extreme size ratios (Table. 2), a high error with the 1-parameter model is expected. This results because of the disparity in surface sites for these ratios relative to the symmetric mixture within which the density independent geometric effects were obtained (Eq. 25). However, even for these extreme size ratios, with information for reference hard sphere distribution functions from simulation, the complete reference theory is able to capture the average bonding numbers for these extreme size ratios quite well.

IV.2.2 Varying association strengths

To understand the effect of varying association strengths between solute-solvent and solvent solvent particles, we studied a case with a fixed size ratio of $\sigma_{s}/\sigma_{p}=0.8$ and for varying association strengths (fig. 7).

As the strength of solute-solvent association is increased as compared to solvent-solvent interactions, multi-body effects become important. It was observed that higher deviations are observed with TPT2-based Marshall-Chapman approximation, especially for increasing strength of solute-solvent association. Importantly, excellent agreement with the complete reference theory is observed for all cases noted in the figure.

IV.2.3 Chemical potential and energy entropy contributions

Fig. 8 shows the chemical potential of the solute for two limiting cases: one where the solvent-solvent association is present and one where it is absent, for two size ratios.

The results indicate that despite the large deviations in $n_{avg}$ (Fig. 6) noted above, the deviations using the Marshall-Chapman approximation are not as high relative to the complete reference approach.

To better understand this result, we consider a symmetric system where there is no association between solvent particles. For this case, the partial molar energy is the energy of the system with solute-solvent association and there is no contribution for change in solvent-solvent association due to inclusion of an ideal solute at infinite dilution. We decompose the residual chemical potential of the solute ($\mu^{Asso}_{s}$) into its energetic and entropic contributions Yu and Karplus (1988); Yu et al. (1990):

[TABLE]

where $E_{s}^{Asso}=(\frac{\partial\mu^{Asso}_{s}/T}{\partial 1/T})_{{\rho_{s}},{\rho_{p}}}$ is the partial molar energy and $S_{s,V}^{Asso}={({\frac{{\partial\mu_{s}^{Asso}}}{{\partial T}}})_{{\rho_{s}},{\rho_{p}}}}$ is the partial molar entropy contribution. The entropy and energy values can be obtained with MCA and complete reference approach based on the corresponding temperature derivatives of the chemical potential of solute.

Fig. 9 clearly shows that energy and entropy contributions are not captured accurately by MCA. But the apparent reasonable prediction of the chemical potential (using MCA) results because of cancellation of the errors between the entropy and energy contributions. The above result shows that in comparing perturbation theories, it could be useful and prudent to study chemical potential and also its energy and entropy contributions.

IV.2.4 Concentrated systems

We also study the variation of average bonding number of the solute for different concentration of the solute in the asymmetric mixtures. As the concentration of the solute increases, the system becomes limited in the number of solvent molecules that can bond to the solute molecules and hence TPT2 approximation used in MCA becomes more accurate. For low concentrations, multi-body correlations are important and hence deviations are observed with MCA. Fig. 10 shows the results for two different size ratios. The theory with the complete reference is able to capture the average bonding numbers for the whole concentration range.

V Concluding discussions

We have studied asymmetric mixtures having strong short-range association between differently sized solute and solvent molecules. The solute molecules have isotropic interactions and the solvent molecules have directional interactions. Such systems are archetypes of colloidal mixtures that are being actively studied in designing materials from the nanoscale level. These systems can also describe the short range ion-solvation and ion-pairing effects in electrolyte systems which is another focus of our research.

The isotropic interactions of the solute can allow multiple solvent molecules to associate and hence multi-body effects become important for these systems. Previously Bansal et al. (2016), we discussed the development of an accurate perturbation theory for these systems is hindered by the difficulty in obtaining the multi-body correlations in the reference system (typically hard sphere). We discussed the limitation of an approach based on obtaining multi-body correlations for the hard sphere system in the gas phase and approximating bulk solvent effect with a linear superposition of pair correlation (together with term to account for three body correction). It was observed that this second order perturbation method fails at high association strengths and high densities. We introduced an approach to represent the multi-body clusters in terms of occupancy distribution to accurately describe the packing in the hard sphere system. Excellent agreement with MC simulation for a range of conditions of association and concentrations were obtained with this complete reference approach for symmetric mixtures.

Here, we have built upon the earlier work and study systems with size asymmetry. Our study shows that the multi-body correlations for asymmetric hard sphere mixtures can be accurately studied in terms of occupancy distributions. With these accurate packing effects, our approach gives excellent agreement with MC simulation for the asymmetric associating mixtures. These occupancy distributions were obtained by particle simulations. Based on ideas borrowed from quasichemical theory, we have also developed parametric models to describe occupancy distributions in the hard sphere systems of different densities and different asymmetries. These distributions were obtained by describing the effects of clustering, medium and surface interactions simultaneously in the hard sphere packing around a solute and can be incorporated in perturbation theories (for eg. statistical associating fluid theory) without having to perform particle simulations for the reference fluid. We validate this complete reference theory (using parameterized models of the reference distribution) against several simulations. A critical test was in analyzing the energy and entropy contributions to the chemical potential of the solute. For a system where solvent-solvent association is not present, such a decomposition of chemical potential in energy and entropy contributions with complete reference theory showed excellent agreement with MC simulations. The apparently reasonable agreement of the second order perturbation approach is shown to arise from the balancing of errors in the energy and entropy contributions. This important finding suggests the need to study different properties while validating perturbative theories for fluids. This present framework can prove useful in modeling real solutions where concentration of solute is low and its size is different from the solvent molecules.

VI Acknowledgment

We acknowledge RPSEA / DOE 10121-4204-01 and the Robert A. Welch Foundation (C-1241) for financial support.

VII Appendix

VII.1 Hard sphere distribution

VII.1.1 Expression for $K_{n}$

From Eq. 19, we find that obtaining an expression for $K_{n}$ reduces to evaluating the ratio of $p_{n}/p_{0}$. The total potential energy of the system when $n$ solvent particles are coordinated with the solute and the remaining $N-n$ solvent particles are outside the observation volume can be formally written as $U=U_{SP_{n}}+U_{N-n|SP_{n}}+U_{N-n}$. $U_{SP_{n}}$ is the potential energy of the solute-$n$-solvent cluster. $U_{N-n|SP_{n}}$ is the interaction energy of the cluster with the rest of the solvent; specifically $U_{N|SP_{0}}$ is the interaction energy of the solute with the fluid outside the observation volume. In the particular case of hard-spheres, $U_{N|SP_{0}}=0$. Finally, $U_{N-n}$ is the potential energy of the solvent constituting the bulk.

Since $p_{n}\propto Q(n,N-n,V,T)$, where $Q(n,N-n,V,T)$ is the canonical partition function of the system with $n$ solvent in the observation volume around the solute and $N-n$ in the bulk, we immediately have

[TABLE]

where we have implicitly moved the center of the coordinates to the center of the solute and thus canceled a common factor of $V$ from both the numerator and denominator. Further, since $U_{SP_{0}}=0$ and $U_{N|SP_{0}}=0$, the denominator simply depends on the potential energy of the solvent in the bulk. (Of course for a general solute, this restriction is easily removed Merchant and Asthagiri (2009); Merchant et al. (2011a).)

Next consider the ratio

[TABLE]

where we suppress $T$ for conciseness and the 0 indicates that there is no solute in the system (or as is the case here, $U_{SP_{0}}=U_{N|SP_{0}}=0$). In the thermodynamic limit of large $V>>v$ and $N>>n$, from the standard potential distribution relation Beck et al. (2006); Pratt and Asthagiri (2007); Widom (1982), each of the above factor on the right is simply $e^{-\beta\mu_{p}^{\rm ex}}/\Lambda_{p}^{3}\rho_{p}$, where $\Lambda_{p}$ is the thermal de Broglie wavelength of the solvent sphere and $\mu^{\rm ex}_{p}$ is its excess chemical potential, and $\rho_{p}$ the density of solvent.

Since,

[TABLE]

we multiply and divide Eq. 28 by the factor $\int\limits_{V-v}d{\vec{r}}_{1}\ldots\int\limits_{V-v}d{\vec{r}}_{N-n}e^{-\beta U_{N-n}}$. Rearranging the resulting equation using Eq. 29 in the large $V$ and large $N$ limit, and noting that the momentum partition functions (for both solute and solvent) cancel exactly, we obtain Eq. 20 (main text), where

[TABLE]

where $\langle\ldots\rangle_{N-n}$ denotes averaging over the configurations of the $N-n$ solvent particles in the volume outside the observation shell.

VII.1.2 MaxEnt model for $\{p_{n}\}$

Here we present an alternative derivation of the two parameter model (Eq. 24) on the basis of information theoretic modeling of $\{p_{n}\}$ Sivia (1996); Hummer et al. (1998); Asthagiri et al. (2003). On the basis of the isolated cluster partition function, we have the distribution of occupancy probabilities $\{p_{n}^{(0)}\}$ as

[TABLE]

With $\{p_{n}^{(0)}\}$ as the default model and accepting the availability of the mean (first moment) and variance (second moment) of the distribution $\{p_{n}\}$ from simulation data, by standard maximum entropy arguments, we have

[TABLE]

where the Lagrange multipliers $C$, $\lambda_{1}$, and $\lambda_{2}$ are, respectively, obtained from enforcing the following constraints

[TABLE]

[TABLE]

[TABLE]

We thus obtain

[TABLE]

Eq. 37 is the same form as obtained in section II.2 for a two parameter correlation. By using $\lambda^{\prime}=e^{-\lambda_{1}}$, we can also represent Eq. 37 by

[TABLE]

which is the same form introduced in our previous work Bansal et al. (2016). Based on this derivation, we can see that $\lambda_{2}$ is the term corresponding to surface interactions discussed in section  II.2 and constraints the variance of the $\{p_{n}\}$ distribution.

VII.2 Isolated cluster probabilities

For asymmetric mixtures we study the isolated cluster probabilities and to find the maximum number of solvent molecules that can occupy the observation volume, we use spherical code Weisstein (1999); Sloane et al. (1994).

VII.2.1 Spherical Code

The spherical code provides information to place $n$ points optimally across a sphere. It gives the optimal angle between the points, considering the center of the sphere as the origin. This is extended to place solvent particles across the solute surface, point of contact between the sphere serving as the optimizing points. The angle between the two contact points, $\theta$, is determined as given in the fig. 11. The number of solvent molecule that can be tightly packed is obtained from the data Weisstein (1999). The sphere onto which the points are optimally placed is an imaginary sphere which includes the critical radius as shown in the dashed lines. It is to be noted that this is still a theoretical estimate for contact packing on the imaginary larger sphere, and due to higher freedom of packing in our case, coordination states can be marginally higher for very large size ratios, $\sigma_{s}/\sigma_{p}\geq 5$. Table in fig. 11 also gives the ma ximum an gle for which single bonding condition holds for a given size ratio ($\theta_{c,max}$) and specified critical distance ($r_{c}$).

Once the $n^{max}$ is defined, we find $P^{(n)}$ as the probability that there is no hard sphere overlap for randomly generated solvent molecules in the observation volume (or inner-shell) of solute molecules. As was discussed previouslyMarshall and Chapman (2013a) that a hit-or-miss Monte Carlo Hammersley and Handscomb (1964); Pratt et al. (2001) approach to calculate $P^{(n)}$ proves inaccurate for large values of $n$ ($n>8$). But since

[TABLE]

where $P_{insert}^{(n)}$ is the probability of inserting a single particle given $n-1$ particles are already in the bonding volume, an iterative procedure was be used to build the higher-order partition function from lower order one Marshall and Chapman (2013a). The one-particle insertion probability $P_{insert}^{(n)}$ is easily evaluated using hit-or-miss Monte Carlo. Following is the table for isolated cluster probabilities for different size ratios.