# Machine-Learning Many-Body Potentials for Colloidal Systems

Gerardo Campos-Villalobos<EMAIL_ADDRESS>Soft Condensed Matter,
Debye Institute for Nanomaterials Science, Utrecht University, Princetonplein
1, 3584 CC Utrecht, The Netherlands  Emanuele Boattini  Soft Condensed Matter,
Debye Institute for Nanomaterials Science, Utrecht University, Princetonplein
1, 3584 CC Utrecht, The Netherlands  Laura Filion  Soft Condensed Matter,
Debye Institute for Nanomaterials Science, Utrecht University, Princetonplein
1, 3584 CC Utrecht, The Netherlands  Marjolein Dijkstra<EMAIL_ADDRESS>Soft
Condensed Matter, Debye Institute for Nanomaterials Science, Utrecht
University, Princetonplein 1, 3584 CC Utrecht, The Netherlands

###### Abstract

Simulations of colloidal suspensions consisting of mesoscopic particles and
smaller species like ions or depletants are computationally challenging as
different length and time scales are involved. Here, we introduce a machine
learning (ML) approach in which the degrees of freedom of the microscopic
species are integrated out and the mesoscopic particles interact with
effective many-body potentials, which we fit as a function of all colloid
coordinates with a set of symmetry functions. We apply this approach to a
colloid-polymer mixture. Remarkably, the ML potentials can be assumed to be
effectively state-independent and can be used in direct-coexistence
simulations. We show that our ML method reduces the computational cost by
several orders of magnitude compared to a numerical evaluation, and accurately
describes the phase behavior and structure, even for state points where the
effective potential is largely determined by many-body contributions.

††preprint: AIP/123-QED

##  I Introduction

Colloidal suspensions consist of mesoscopic particles suspended in a solvent.
The effective interactions between the colloids can be modified by the
addition of smaller species like ions, depletants, ligands, or polymers. For
instance, the addition of non-adsorbing polymer to a colloidal suspension
induces an effective attraction between the colloids. This arises from an
increase in configurational entropy of the polymer as colloidal particles
approach each other and their depletion zones overlap as described by Asakura
and Oosawa in 1954. Asakura and Oosawa (1954) The strength and range of this
depletion interaction can be independently tuned by varying the polymer
fugacity and polymer size. Gast, Hall, and Russel (1983); Meijer and Frenkel
(1991); Lekkerkerker  _et al._ (1992); Ilett  _et al._ (1995); Dijkstra,
Brader, and Evans (1999) The possibility of tailoring the effective
interactions enriches the physics of colloidal systems compared to simple
atomic fluids, and leads to a wide range of potential applications.

Figure 1:  Left: Schematic configuration of a colloid (blue sphere) with its
corresponding depletion zone (gray-shaded sphere). Right: The values of the
descriptor $n=n(\bm{r})$, indicating the number of simultaneously overlapping
depletion layers, are drawn for different regions in the triplet
configuration. The region highlighted in pink, with $n=3$, corresponds to the
three-body correction, $V_{\text{f}}^{(3+)}$, that is added on top of the
pair-wise approximation of the total free volume. In general, non-zero
contributions to $V_{\text{f}}^{(3+)}$ will stem from regions where $n\geq 3$.

Computer simulations of these systems remain computationally expensive due to
slow dynamics as very different length and time scales are involved for the
various species and as the number of microscopic species outweighs by orders
of magnitude the number of colloidal particles. Dijkstra (2001) Huge efforts
have been devoted to speeding up equilibration in highly asymmetric mixtures
by the implementation of (rejection-free) cluster moves, Dijkstra and van Roij
(1997); Buhot and Krauth (1998); Lobaskin and Linse (1999); Liu and Luijten
(2004); Vink and Horbach (2004); Bernard, Krauth, and Wilson (2009) lattice
discretization methods, Meijer and Frenkel (1991); Dijkstra and Frenkel
(1994); Panagiotopoulos and Kumar (1999) or by exploiting a sophisticated
field-theoretic description of the smaller species in the external field of a
fixed colloid configuration within a density functional framework Löwen,
Madden, and Hansen (1992); Fushiki (1992) in the same spirit as the "ab
initio" method of Car and Parrinello for ion-electron systems. Car and
Parrinello (1985) However, all these approaches have their own limitations,
and as a consequence higher densities, larger system sizes and larger size
asymmetries are still unattainable, thereby leaving many intriguing
experimental observations unexplained such as void formation, gas-liquid and
gas-solid coexistence in like-charged colloidal suspensions, Tata,
Rajalakshmi, and Arora (1992); Ito, Yoshida, and Ise (1994); Kepler and Fraden
(1994); Tata  _et al._ (1997); Larsen and Grier (1997) as well as hampering
investigations of interesting phenomena like capillary wave fluctuations, long
wavelength fluctuations of marginal colloidal liquids at their triple point,
and density oscillations at the gas-liquid interface in colloid-polymer
mixtures. Aarts, Schmidt, and Lekkerkerker (2004); Vink, Horbach, and Binder
(2005); Moussaid  _et al._ (1999); Vis  _et al._ (2020)

Figure 2:  Square of the correlation coefficient $R^{2}$ and root mean squared
error RMSE as a function of the number of selected SFs $N_{\text{s}}$ for a
colloid-polymer mixture with size ratio $q=0.4$, $0.8$ and $1.0$. The RMSE is
shown for both the training and test sets. The optimal value of the number of
SFs is indicated with the vertical dashed line for each case. Note that
$R^{2}$ is plotted on linear scale and RMSE on logarithmic scale.

An alternative strategy to circumvent slow equilibration is to formally
integrate out the degrees of freedom of the microscopic species in the
partition sum and to derive an exact expression for the effective one-
component Hamiltonian of the colloids. The effective Hamiltonian, which
depends on all colloid coordinates and involves many-body interactions, can be
employed in standard simulation schemes, but its evaluation becomes extremely
computationally demanding when three- and higher-body interactions are
important. Dijkstra  _et al._ (2006); Kobayashi  _et al._ (2019)

In the past years, machine learning (ML) techniques have been employed to
efficiently approximate the many-body interatomic potentials in atomistic
systems by fitting large data sets from electronic structure calculations.
Rupp  _et al._ (2012); Behler (2016); Musil  _et al._ (2021) These ML
potentials can be used in molecular simulations, thereby combining the
accuracy of “first principles" calculations with the efficiency of simple
atomistic “force-fields". Behler (2016)

In this Article, we introduce a general ML approach to efficiently represent
the effective many-body interactions in colloidal systems. We consider a model
suspension of sterically-stabilized colloidal particles and non-adsorbing
ideal polymers for which a wealth of data exist on the phase behavior and
structure, both from numerical evaluations of the exact one-component
Hamiltonian or from direct simulations of the true binary mixture.Dijkstra,
Brader, and Evans (1999); Dijkstra  _et al._ (2006); Verso  _et al._ (2006);
Vink and Horbach (2004) Moreover, exemplifying our method using this model
system is convenient as the importance of two-, three-, and higher-body
contributions to the effective potential can be tuned by the polymer size and
colloid density, allowing for a systematic exploration of the validity of our
ML approach. We fit the effective many-body potentials as a function of all
colloid coordinates with the symmetry functions (SFs) as introduced by Behler
and Parrinello, Behler and Parrinello (2007) which were also recently used to
fit the many-body interactions between microgel particles. Boattini  _et al._
(2020) We find that the fitted ML potential reduces the computational cost by
up to almost four orders of magnitude in comparison to a numerical evaluation,
and accurately describes the effective many-body interactions for a wide range
of colloid densities and polymer fugacities. Hence, the ML potential can be
assumed to be state-independent and can straightforwardly be used in Monte
Carlo (MC) simulations to observe direct coexistence of a dilute colloidal gas
and a dense colloidal liquid phase. Using these ML potentials, we find good
agreement with the structure, phase behavior, and interfacial tension as
obtained in previous studies for mixtures of different size ratio.

The remainder of this article is structured as follows. In Sec. II, we
describe the model colloid-polymer mixture and its formal mapping onto the
corresponding effective one-component representation. The employed SFs and
fitting procedure of the ML potentials are presented in Sec. III, whereas the
results of the MC simulations using the ML potentials, demonstrating the
accuracy of the model, are discussed in Sec. IV. We conclude with a final
discussion in Sec. V.

##  II Effective one-component description of the Asakura-Oosawa model

The simplest model that captures the essence of polymer-induced effective
interactions between colloidal particles Asakura and Oosawa (1954, 1958) was
introduced by Vrij. Vrij (1976) In this so-called Asakura-Oosawa (AO) model,
the colloidal particles are regarded as hard spheres of diameter
$\sigma_{\text{c}}$, whereas the polymer coils with radius of gyration
$R_{\text{g}}$ are treated as ideal point particles as regards their mutual
interactions. The colloid-polymer pair interaction is hard-sphere-like such
that their distance of closest approach is
$\sigma_{\text{cp}}=(\sigma_{\text{c}}+\sigma_{\text{p}})/2$ with polymer
diameter $\sigma_{\text{p}}=2R_{\text{g}}$. We consider a system consisting of
$N_{\text{c}}$ colloidal hard spheres at coordinates $\\{\bm{R}_{i}\\}$ with
$i=1,\dots,N_{\text{c}}$, and $N_{\text{p}}$ polymer coils at positions
$\\{\bm{r}_{j}\\}$ with $j=1,\dots,N_{\text{p}}$ with a size ratio
$q=\sigma_{\text{p}}/\sigma_{\text{c}}$ in a volume $V$ at temperature $T$.
The pair interactions of the AO model read

$\phi_{\text{cc}}(R_{ij})=\begin{cases}\infty&\text{for
}R_{ij}<\sigma_{\text{c}}\\\ 0&\text{otherwise, }\end{cases}$ (1)
$\phi_{\text{cp}}(|\bm{R}_{i}-\bm{r}_{j}|)=\begin{cases}\infty&\text{for
}|\bm{R}_{i}-\bm{r}_{j}|<\sigma_{\text{cp}}\\\ 0&\text{otherwise,
}\end{cases}$ (2) $\phi_{\text{pp}}(r_{ij})=0,$ (3)

with $R_{ij}=|\bm{R}_{i}-\bm{R}_{j}|$ the center-of-mass distance between
colloid $i$ and $j$. This binary mixture of colloids and polymer is described
by the total interaction Hamiltonian
$\mathcal{H}=\mathcal{H}_{\text{cc}}+\mathcal{H}_{\text{cp}}+\mathcal{H}_{\text{pp}}$
with

$\mathcal{H}_{\text{cc}}=\sum_{i<j}^{N_{\text{c}}}\phi_{\text{cc}}(R_{ij}),$
(4)
$\mathcal{H}_{\text{cp}}=\sum_{i=1}^{N_{\text{c}}}\sum_{j=1}^{N_{\text{p}}}\phi_{\text{cp}}(|\bm{R}_{i}-\bm{r}_{j}|),$
(5) $\mathcal{H}_{\text{pp}}\equiv 0.$ (6)

Figure 3:  The probability of $n$ overlapping depletion layers for a colloid-
polymer mixture with polymer reservoir packing fraction
$\eta_{\text{p}}^{\text{r}}=0$ and colloid packing fraction $\eta_{\text{c}}$
as labelled for varying size ratios $q=0.4$, $0.8$ and $1.0$.

It is convenient to treat the polymer coils grand-canonically, in which the
fugacity of the polymers $z_{\text{p}}$, or equivalently the polymer reservoir
packing fraction
$\eta_{\text{p}}^{\text{r}}\equiv\pi\sigma_{\text{p}}^{3}z_{\text{p}}/6$, is
fixed. The thermodynamic potential $F(N_{\text{c}},z_{\text{p}},V,T)$ of this
binary system reads

$\displaystyle\exp[-\beta F]$ $\displaystyle=$
$\displaystyle\sum_{N_{\text{p}}=0}^{\infty}\frac{z_{\text{p}}^{N_{\text{p}}}}{N_{\text{c}}!\Lambda_{\text{c}}^{3N_{\text{c}}}N_{\text{p}}!}\int_{V}d\bm{R}^{N_{\text{c}}}\int_{V}d\bm{r}^{N_{\text{p}}}\exp[-\beta(\mathcal{H}_{\text{cc}}+\mathcal{H}_{\text{cp}})],$
(7) $\displaystyle=$
$\displaystyle\frac{1}{N_{\text{c}}!\Lambda_{\text{c}}^{3N_{\text{c}}}}\int_{V}d\bm{R}^{N_{\text{c}}}\exp\left[-\beta\mathcal{H}_{\text{cc}}\right]\exp\left[-\beta\Omega\right]$
$\displaystyle=$
$\displaystyle\frac{1}{N_{\text{c}}!\Lambda_{\text{c}}^{3N_{\text{c}}}}\int_{V}d\bm{R}^{N_{\text{c}}}\exp[-\beta\mathcal{H}_{\text{eff}}],$

where $\beta=(k_{\text{B}}T)^{-1}$ with $k_{\text{B}}$ the Boltzmann constant,
and $\Lambda_{\alpha}$ is the thermal wavelength of species
$\alpha=\text{c},\text{p}$. By integrating out the degrees of freedom of the
polymers, this binary mixture can be mapped onto an effective one-component
system described by an effective colloids-only Hamiltonian
$\mathcal{H}_{\text{eff}}=\mathcal{H}_{\text{cc}}+\Omega$, where
$\Omega=-\beta^{-1}\ln\left[\sum_{N_{\text{p}}=0}^{\infty}\frac{z_{\text{p}}^{N_{\text{p}}}}{N_{\text{p}}!}\int_{V}d\bm{r}^{N_{\text{p}}}\exp\left(-\beta\mathcal{H}_{\text{cp}}\right)\right]$
denotes the grand potential of a “sea" of ideal polymers at fugacity
$z_{\text{p}}$ in the external field of a fixed configuration of
$N_{\text{c}}$ colloids. Dijkstra, Brader, and Evans (1999); Dijkstra  _et
al._ (2006); Dijkstra and van Roij (2002) For the AO model, the grand
potential $\Omega$ is simply the negative of the free volume available for the
polymer in the fixed configuration of $N_{\text{c}}$ colloids, i.e.

$\Omega=-z_{\text{p}}V_{\text{f}}(\\{\bm{R}_{i}\\}),$ (8)

which is the volume outside the $N_{\text{c}}$ depletion zones, and which can
be decomposed into a zero-, one-, two-, three-, and higher-body contribution

$V_{\text{f}}=V_{\text{f}}^{(0)}+\sum_{i=1}^{N_{\text{c}}}{V_{\text{f}}^{(1)}(\bm{R}_{i})}+\sum_{i<j}^{N_{\text{c}}}{V_{\text{f}}^{(2)}(\bm{R}_{i},\bm{R}_{j})}+V_{\text{f}}^{(3+)},$
(9)

i.e. the number of colloids $k=0,1,2,\cdots,N_{\text{c}}$ that interact
simultaneously with the “sea” of ideal polymer, and where
$V_{\text{f}}^{(3+)}$ denotes the three- and higher-body term. Analytical
expressions exist for $k=0,1,$ and 2. The zero-body contribution
$V_{\text{f}}^{(0)}=V$ simply corresponds to the volume of the system,
$V_{\text{f}}^{(1)}=-v_{1}$ is the volume excluded to a polymer by a single
colloid with $v_{1}=\pi\sigma_{\text{cp}}^{3}/6$, and
$V_{\text{f}}^{(2)}=8v_{1}(1-3x/4+x^{3}/16)$ with
$x=R_{ij}/\sigma_{\text{cp}}<2$ is the well-known depletion potential of the
AO model, representing the lens-shaped overlap volume of two colloidal spheres
of diameter $\sigma_{\text{cp}}$ at separation $R_{ij}$.

For $q<0.1547$, the higher-order contributions incorporated in
$V_{\text{f}}^{(3+)}$ are zero, and hence a mapping onto an effective one-
component system with an effective Hamiltonian based on pairwise additive
depletion potentials is exact. Gast, Hall, and Russel (1983); Dijkstra,
Brader, and Evans (1999) For larger $q$, the model based solely on the
pairwise approximation exhibits a phase behaviour that strongly differs from
that observed in the system where three- and higher-body interactions
$V_{\text{f}}^{(3+)}$ are considered. The main effect of the many-body
interactions in colloid-polymer mixtures is to enhance the $z_{\text{p}}$
regime of stable gas-liquid coexistence at the expense of that of the gas-
solid coexistence. Dijkstra, Brader, and Evans (1999); Dijkstra  _et al._
(2006) $V_{\text{f}}^{(3+)}$ can be computed by measuring the overlap volume
of three or more depletion zones which can only be evaluated numerically as
shown in Ref. Dijkstra  _et al._ , 2006. The key ingredient in the numerical
evaluation of $V_{\text{f}}^{(3+)}$ in Ref. Dijkstra  _et al._ , 2006 is the
introduction of a spatial descriptor $n=n(\bm{r})$, which for a given colloid
configuration, measures the number of simultaneously overlapping depletion
layers at spatial coordinate $\bm{r}$. In Fig. 1 we show a pictorial
representation of a configuration of colloidal particles (solid blue spheres)
and their depletion zones (surrounding gray-shaded spheres). The values of $n$
in different regions of space are also reported. Upon the introduction of
$n(\bm{r})$, the many-body correction $V_{\text{f}}^{(3+)}$ can be computed as

$V_{\text{f}}^{(3+)}=-\frac{1}{2}\int_{n\geq
3}{d\bm{r}\left[n(\bm{r})-1\right]\left[n(\bm{r})-2\right]},$ (10)

where the integration is performed over the regions with $n\geq 3$. Dijkstra
_et al._ (2006) Conceptually, $V_{\text{f}}^{(3+)}$ can be seen as a
correction to the overestimated overlap volume of the depletion zones in
$V_{\text{f}}^{(2)}$.

##  III Fitting the many-body potential

As described above, $V_{\text{f}}^{(3+)}$ is a function of the colloid
configuration $\bm{R}^{N_{\text{c}}}$. Therefore, its evaluation requires to
sample the whole volume $V$ in order to identify those regions where $n\geq
3$. In practice, $n(\bm{r})$ can be measured using a spherically-symmetric
grid of $M$ points around each colloid. The number of points needs to be
sufficiently high (on the order of $M\sim 10^{5}$) to accurately sample
$V_{\text{f}}^{(3+)}$, making its implementation computationally feasible only
for a few hundreds of colloids. Here, we use an ML approach to fit
$V_{\text{f}}^{(3+)}$ as a function of all colloid coordinates using a set of
SFs introduced by Behler and Parrinello. Behler (2011); Behler and Parrinello
(2007)

###  III.1 Training data set

Since $V_{\text{f}}^{(3+)}$ depends on size ratio $q$, we use a different ML
fit for each $q$. To build the training data sets, we perform MC simulations
on $N_{\text{c}}=108$ colloids for $q=0.4,0.8$, and 1.0, polymer fugacity
$z_{\text{p}}=0$, and colloid packing fraction
$\eta_{\text{c}}=\pi\sigma_{\text{c}}^{3}N_{\text{c}}/6V\in\left[0.15,0.65\right]$
with a packing fraction spacing of $\delta\eta_{\text{c}}=0.005$. From each
simulation we collect 500 equilibrated, essentially uncorrelated
configurations111The configurations we use are spaced far enough apart in
number of MC cycles that the particles have had a chance to diffuse
significantly between successive snapshots. Even in the high density fluid
phase at $\eta_{\text{c}}=0.40$, the particles travel more than 6 times their
diameters between successive configurations. and measure $V_{\text{f}}^{(3+)}$
using a spherically-symmetric $(r^{3},\cos(\theta),\phi)$ grid of 100, 50 and
50 points, respectively. Dijkstra  _et al._ (2006) The resulting data set for
each $q$ contains a total of $49,500$ representative particle configurations
at different colloid densities, from which $80\%$ are used for training and
$20\%$ for testing. It is worth mentioning that the selected case of zero
polymer fugacity corresponds to a single-component system of hard spheres,
which is computationally extremely inexpensive to simulate. Therefore,
equilibration in this system is achieved rapidly and we were able to access a
large number of decorrelated configurations with relevant local environments
of the individual particles in low- and high-density fluid and solid phases,
which are also present in colloid-polymer mixtures. Learning
$V_{\text{f}}^{(3+)}$ instead of the full grand potential $\Omega$
significantly reduces the size of the training data set and lowers the
computational effort while still giving access to non-zero polymer fugacity
via Eq. 8.

Figure 4:  Ratio between the computational time required for the numerical
evaluation of the many-body term term ($\tau_{\text{N}}$) and the ML potential
($\tau_{\text{ML}}$) as a function of colloid packing fraction
$\eta_{\text{c}}$ of a colloid-polymer mixture with varying size ratio $q$ and
polymer fugacity $z_{\text{p}}=0$. Figure 5:  The two-body
$V_{\text{f}}^{(2)}$ (pentagons), and three- and higher-body
$V_{\text{f}}^{(3+)}$ contribution to the effective many-body potential as
predicted by the ML potential (ML) (squares) for a colloid-polymer mixture
with size ratio $q=0.4,0.8$, and 1.0 as a function of colloid packing fraction
$\eta_{\text{c}}$ and varying polymer reservoir packing fraction
$\eta_{\text{p}}^{\text{r}}$ along with the actual $V_{\text{f}}^{(3+)}$ as
obtained from a numerical evaluation (N) (plus). Inset plots correspond to the
percentage errors of the many-body term, $\Delta
V_{\text{f}}^{(3+)}=\left[|V_{\text{f}}^{(3+)}(\text{N})-V_{\text{f}}^{(3+)}(\text{ML})|\right]/V_{\text{f}}^{(3+)}(\text{N})\times
100$. Note that for $\eta_{\text{p}}^{\text{r}}\neq 0$ we consider only state
points outside coexistence regions. Figure 6:  Colloid-colloid (squares),
colloid-polymer (pentagons) and polymer-polymer (circles) pair correlation
functions $g_{\alpha\beta}(r)$ with $\alpha,\beta=\text{c},\text{p}$ of a
colloid-polymer mixture with $q=0.4$, $0.8$ and $1.0$ obtained from MC
simulations using the ML potential at $\eta_{\text{c}}$ and
$\eta_{\text{p}}^{\text{r}}$ as labeled. In the top and middle panels, the
solid lines are the results obtained from MC simulations based on a numerical
evaluation of $V_{\text{f}}^{(3+)}$, whereas in the bottom panel, the solid
lines correspond to the results reported in Ref. Dijkstra  _et al._ , 2006.
For clarity, $g_{\text{cp}}(r)$ and $g_{\text{cc}}(r)$ are shifted in the
vertical direction.

###  III.2 Symmetry functions

To describe the local environment of a particle we use the SFs introduced by
Behler and Parrinello for constructing high-dimensional neural network
potentials. Behler and Parrinello (2007) These SFs are described in great
detail in Refs. Behler, 2011, 2015 and have been used as inputs for atomic
feed-forward neural networks in order to provide the atomic energy
contributions of different materials and molecules. Behler and Parrinello
(2007); Khaliullin  _et al._ (2010); Eshet  _et al._ (2012); Kapil, Behler,
and Ceriotti (2016); Cheng, Behler, and Ceriotti (2016) Recently, they have
also been used in combination with linear regression to fit effective many-
body interactions between elastic spheres. Boattini  _et al._ (2020)

Since the three- and higher-body term $V_{\text{f}}^{(3+)}$ in Eq. 10 does not
include any two-body contributions, we consider only the angular three-body
SFs $G^{3}(i)$ for particle $i$, which are defined as

$\displaystyle G^{3}(i)=$ $\displaystyle 2^{1-\xi}\sum_{j,k\neq
i}\left(1+\lambda\cos\theta_{ijk}\right)^{\xi}e^{-\eta\left(R_{ij}^{2}+R_{ik}^{2}+R_{jk}^{2}\right)}$
(11) $\displaystyle\times f_{c}(R_{ij})f_{c}(R_{ik})f_{c}(R_{jk}),$

where the indices $j$ and $k$ run over all the neighbours of particle $i$, and
$\xi$, $\eta$, and $\lambda$ are three parameters that determine the shape of
the function. The parameter $\lambda$ can have the values $+1$ or $-1$ and
determines the angle $\theta_{ijk}$ at which the angular part of the function
has its maximum. The angular resolution is provided by the parameter $\xi$,
while $\eta$ controls the radial resolution. Additionally, $f_{c}(R_{ij})$ is
a cutoff function: a monotonically decreasing function that smoothly goes to 0
in both value and slope at the cutoff distance $r_{c}$. Here, we consider a
cutoff function of the form

$f_{c}(R_{ij})=\begin{cases}\tanh^{3}(1-R_{ij}/r_{c})&\text{for }R_{ij}\leq
r_{c}\\\ 0&\text{for }R_{ij}>r_{c}.\end{cases}$ (12)

Using the training set consisting of representative colloid configurations at
different colloid densities, we fit $V_{\text{f}}^{(3+)}$ with a linear
combination of $N_{\text{s}}$ SFs in Eq. 11, and select the optimal subset of
SFs using the feature selection scheme of Ref. Boattini  _et al._ , 2020,
which we summarize in the following. We note that we also implemented other
feature selection schemes such as those introduced by Imbalzano et al.,
Imbalzano  _et al._ (2018) but the one of Ref. Boattini  _et al._ , 2020 turns
out to be more efficient (see Appendix).

The first step of the method involves the creation of a large but manageable
pool of candidate SFs. This is done by calculating, for every configuration in
the data set, several SFs with different sets of parameters. Specifically, we
generate the $G^{3}(i)$ SFs by setting $r_{c}=2\sigma_{\text{cp}}$,
$\lambda\in\\{-1,1\\}$, $\eta\in\\{0.001,0.01,0.1,1,2,4,8\\}$, and
$\xi\in\\{1,2,4,8\\}$. With these choices, our pool of candidates consists of
$D=56$ SFs.

Then, an optimal subset of $N_{\text{s}}<D$ SFs is selected from the pool in a
step-wise fashion. First, the SF is selected with the largest correlation with
the many-body term as quantified by the Pearson correlation coefficient
$c_{k}$

$c_{k}=\frac{\sum_{j}{\left(\sum_{i}G_{k}^{3}(i)|_{j}-\overline{\sum_{i}G_{k}^{3}(i)}\right)\left(V_{\text{f}}^{(3+)}|_{j}-\overline{V_{\text{f}}^{(3+)}}\right)}}{\sigma_{\text{SD}}(\sum_{i}G_{k}^{3}(i))\sigma_{\text{SD}}(V_{\text{f}}^{(3+)})},$
(13)

where $\sum_{i}G_{k}^{3}(i)|_{j}$ represents the sum of the $k$-th SF over all
colloidal particles $i$ in configuration $j$ and $V_{\text{f}}^{(3+)}|_{j}$
denotes the many-body correction evaluated for this configuration.
$\overline{\sum_{i}G_{k}^{3}(i)}$ and $\overline{V_{\text{f}}^{(3+)}}$
correspond to arithmetic means over the whole data set, and
$\sigma_{\text{SD}}(\sum_{i}G_{k}^{3}(i))$ and
$\sigma_{\text{SD}}(V_{\text{f}}^{(3+)})$ to their standard deviations. The
next SF is then selected based on the highest increase in the linear
correlation between the currently selected set and the target many-body term
as determined by

$R^{2}={\bf c}^{T}{\bf R}^{-1}{\bf c},$ (14)

where ${\bf c}^{T}=(c_{1},c_{2},\cdots)$ is the vector whose $j$-th component
is given by the Pearson correlation coefficient (13) between the $j$-th SF and
the many-body term, and ${\bf R}$ is the correlation matrix of the current set
of SFs with elements $\mathcal{R}_{ij}$ representing the Pearson correlation
function between the $i$-th and $j$-th SF. This choice guarantees that only
SFs that add relevant information are selected, while penalizing highly
correlated SFs with only redundant information as well as SFs that are
sensitive to aspects of the particle’s environment that poorly correlate with
the target many-body term. This process is repeated iteratively and new SFs
are selected until the correlation stops increasing appreciably. Finally, the
selected subset of SFs is used to approximate the target many-body term via
simple linear regression.

###  III.3 Accuracy of the ML potentials

In Fig. 2 we report the correlation coefficient $R^{2}$ and the root mean
squared error (RMSE) of the linear fits with the actual $V_{\text{f}}^{(3+)}$
as a function of the number of SFs for $q=0.4,0.8,$ and 1.0, for both the
training and the test set. We note that RMSE and $R^{2}$ are related by a
simple relation
$R^{2}=1-\text{RMSE}^{2}/\sigma^{2}_{\text{SD}}(V_{\text{f}}^{(3+)})$. Upon
increasing $q$, we clearly observe that an increasing number of SFs is
required to accurately approximate $V_{\text{f}}^{(3+)}$, which can be
understood as the thickness of the depletion layers increases with $q$,
thereby enhancing the many-body effects. To quantify the importance of the
many-body contributions to the effective potential, we calculate $P(n)$, the
probability that we find $n=n(\bm{r})$ overlapping depletion layers at spatial
coordinate $\bm{r}$ in a system of $N_{\text{c}}$ colloids in a volume $V$ at
polymer fugacity $z_{\text{p}}$, and size ratio $q$. In Fig. 3 we show $P(n)$
for the three considered size ratios $q=0.4,0.8$, and 1.0 at varying packing
fraction $\eta_{\text{c}}$ and $\eta_{\text{p}}^{\text{r}}=0$. For $q=0.4$, we
find that the largest number $n$ with non-zero probability is $n=4$, and hence
the effective potential consists of up to four-body contributions. For $q=0.8$
and $1.0$, we find that 6- and 7-body contributions become non-negligible at
high colloid densities. We choose $N_{\text{s}}=12,23,$ and 25 for
$q=0.4,0.8,$ and 1.0, respectively, for which we find good agreement between
the ML fits and the actual many-body $V_{\text{f}}^{(3+)}$ term.

###  III.4 Efficiency of the ML potentials

To quantify the efficiency of our ML potentials, we determine the ratio
between the computational time for the numerical evaluation of the many-body
term ($\tau_{\text{N}}$) and the ML potential ($\tau_{\text{ML}}$) for
different values $\eta_{\text{c}}$ and $q$. The ratios
$\tau_{\text{N}}/\tau_{\text{ML}}$ are extracted by selecting 10 decorrelated
configurations for each colloid packing fraction $\eta_{\text{c}}$ from the
training data set, evaluating the serial computing times for the calculation
of the many-body term and repeating this procedure 100 times to get the
average ratios. Note that the implementations of both algorithms was simple
and no neighbour lists were used. Moreover, both codes were compiled with the
same compiler optimizations. From Fig. 4, we find that our ML potentials speed
up the many-body term evaluation at least by two orders of magnitude for size
ratio $q=1.0$ up to almost four orders of magnitude for $q=0.4$. The speed up
decreases with colloid density $\eta_{\text{c}}$ (due to the increasing number
of neighbours) and with polymer size as the cut-off value of the interaction
depends on $q$. The observed speed up achieved with our ML potentials can be
rationalized by comparing the order of the computations on which both
algorithms are based. When the many-body correction is evaluated in a system
of $N_{\text{c}}$ colloids through numerical integration, the time needed for
such a computation scales as $N_{\text{c}}^{2}M$, with $M$ the number of grid
points. In contrast, if no tricks are used, the time required for the
evaluation of $V_{\text{f}}^{(3+)}$ using the three-body SFs scales as
$N_{\text{c}}^{3}$. Thus, as it is verified in our case, a significant speed
up using the ML potentials can be achieved as long as $N_{\text{c}}<M$.

Figure 7:  Typical configuration of a colloidal gas-liquid phase coexistence
of a colloid-polymer mixture as obtained from direct coexistence simulations
using the ML potential. The blue solid spheres represent the colloidal
particles, whereas the gray layers surrounding the colloids illustrate their
depletion zones. The polymers are integrated out and invisible. Figure 8:
Colloidal gas-liquid binodals in the colloid packing fraction
$\eta_{\text{c}}$ \- polymer reservoir packing fraction
$\eta_{\text{p}}^{\text{r}}$ plane of a colloid-polymer mixture with size
ratio $q=0.4$, $0.8$ and $1.0$ as obtained from direct-coexistence simulations
using the ML potentials (open squares) and from earlier reports (solid
circles), either by Dijkstra et al. Dijkstra  _et al._ (2006) using a
numerical evaluation of the effective many-body interactions or by Lo Verso et
al. Verso  _et al._ (2006) employing simulations of the true binary mixture.
The dashed lines represent the binodals predicted by free-volume theory
(FVT).Dijkstra, Brader, and Evans (1999); Dijkstra  _et al._ (2006) Inset
plots correspond to the percentage error of the average density of the
coexisting gas (solid pentagons) and liquid (empty circles) phases at fixed
$\eta_{\text{p}}^{\text{r}}$,
$\Delta\eta_{\text{c}}=\left[|\eta_{\text{c}}(\text{N})-\eta_{\text{c}}(\text{ML})|\right]/\eta_{\text{c}}(\text{N})\times
100$.

##  IV Validation

In the following, we focus on the validation of the effective one-component ML
many-body potential. It is important to note that we fit the
$V_{\text{f}}^{(3+)}$ term for a fixed $q$ with a single ML potential for all
colloid packing fractions $\eta_{\text{c}}\in\left[0.15,0.65\right]$ and
$\eta_{\text{p}}^{\text{r}}=0$. We now investigate the transferability of the
ML potentials to state points outside the training set by considering colloid
configurations at finite non-zero polymer fugacities $z_{\text{p}}$. In Fig.
5, we show the ensemble average of $V_{\text{f}}^{(3+)}$ of a colloid-polymer
mixture with size ratio $q=0.4,0.8$, and 1.0 as a function of colloid packing
fraction $\eta_{\text{c}}$ and varying polymer reservoir packing fraction
$\eta_{\text{p}}^{\text{r}}$ as obtained from independent MC simulations of
$N_{\text{c}}=108$ colloidal particles using the ML potential along with the
actual $V_{\text{f}}^{(3+)}$ as obtained from a numerical evaluation. The
values of $V_{\text{f}}^{(2)}$ are also included to appreciate the importance
(magnitude) of the many-body and two-body terms at different state points. We
find good agreement between the ML predictions and the actual many-body
$V_{\text{f}}^{(3+)}$ term for the whole range of colloid packing fractions
$\eta_{\text{c}}$ and polymer fugacities $z_{\text{p}}$, demonstrating the
transferability to finite $z_{\text{p}}$ outside the training set. Some
deviations are observed for $q=1.0$ at high $\eta_{\text{c}}$, where four- and
higher-body interactions are predominant. Hence, the ML potential can be
assumed to be effectively state-independent, allowing us to simulate direct
coexistence of an extremely dilute colloidal gas phase with a very dense
colloidal liquid phase in a single simulation box, see Fig. 7.

Figure 9:  Colloid packing fraction profiles for the free interface between
coexisting gas and liquid phases of a colloid-polymer mixture with size ratio
$q=0.8$ at different polymer reservoir packing fraction
$\eta_{\text{p}}^{\text{r}}$, showing weak oscillations at the "colloid-rich"
liquid side near the triple point (a). Gas-liquid interfacial tension
$\beta\gamma\sigma_{\text{c}}^{2}$ for the same mixture as a function of the
difference in gas and liquid packing fractions
$\eta_{\text{c}}^{\text{L}}-\eta_{\text{c}}^{\text{G}}$ as obtained from
direct-coexistence simulations using the ML potential (open squares) and
reported by Vink et al. Vink and Horbach (2004) (solid circles) (b).

Given that the ML potentials accurately describe the many-body
$V_{\text{f}}^{(3+)}$ term at different state points, we now test the ability
of the ML potential in reproducing the structure of a colloid-polymer mixture.
To this end, we measure the colloid-colloid, polymer-polymer and colloid-
polymer pair-correlation functions $g_{\alpha\beta}(r)$ with
$\alpha,\beta=\text{c},\text{p}$ from MC simulations of $N_{\text{c}}=1372$
particles interacting with the ML potential at three different state points
for each size ratio $q=0.4$, 0.8 and 1.0. At first sight it may seem
surprising that it is possible to recover information about the structure of
the polymer as we traced out the polymer degrees of freedom. However, as the
polymer coils are ideal, the number density of the polymer is constant in the
free volume (or holes) of the system.Dijkstra  _et al._ (2006) Therefore, we
can determine the colloid-polymer and polymer-polymer pair correlation
functions from $10,000$ randomly inserted polymer coils in an instantaneous
colloid configuration, provided that no overlap exists between the polymer and
colloids. In Fig. 6 we plot $g_{\alpha\beta}(r)$ denoted by the symbols for
the three size ratios considered at varying state points along with the ones
obtained by a numerical evaluation of $V_{\text{f}}^{(3+)}$ or from previous
work Dijkstra  _et al._ (2006) as denoted by the solid lines. The agreement
between the pair correlation function as obtained using the ML potentials with
the "exact" results is evident in all cases, especially for mixtures with
$q=0.4$ and 0.8, which reflects the ability of our ML potential to correctly
describe the structure of the mixture, regardless of
$\eta_{\text{p}}^{\text{r}}$ and $\eta_{\text{c}}$. However, it is important
to note that mapping of four- and higher-body potentials onto a linear
combination of angular three-body SFs does not necessarily guarantee that the
structure of the colloid-polymer mixture is well captured, see e.g. Refs.
Dijkstra, van Roij, and Evans, 2000; Louis, 2002. Indeed, small deviations are
observed for the contact value of the colloid-colloid pair correlation
function $g_{cc}(r)$ for the highest $q=1.0$, where four- and higher-body
interactions become dominant.

As the ML potentials accurately describe $V_{\text{f}}^{(3+)}$ and the
structure of colloid-polymer mixtures, we now investigate whether the ML
potentials can be used to determine the coexistence densities of the colloidal
gas and liquid phases using direct-coexistence simulations. To this end, we
perform MC simulations of $N_{\text{c}}=1372$ particles in an elongated box.
We determine the packing fraction of the coexisting gas and liquid phases,
$\eta_{\text{c}}^{\text{L}}$ and $\eta_{\text{c}}^{\text{G}}$, by measuring
the equilibrium density profiles. We plot the resulting gas-liquid binodals in
Fig. 8 for colloid-polymer mixtures with size ratio $q=0.4,0.8$, and 1.0. For
$q=0.4$ and $0.8$, we compare our results with those obtained from grand-
canonical MC simulations using the full binary mixture, Verso  _et al._ (2006)
whereas the results for $q=1.0$ are compared with those obtained using the
full effective one-component system. Dijkstra  _et al._ (2006) Our ML results
show good agreement with the previously obtained gas-liquid binodals from the
critical point at low $\eta_{\text{p}}^{\text{r}}$ all the way to the triple
point at high $\eta_{\text{p}}^{\text{r}}$, even for $q=1.0$ where the four-
and higher-body terms contribute for more than 75% to the effective potential
at the triple point as shown in Fig. 3. This is an interesting result, as we
expect a reduced correspondence at large $q$ as the number of depletion layers
that can simultaneously overlap increases due to larger depletion layers. In
addition, we plot the results from free-volume theory for $q=0.8$ and 1.0 and
find good agreement, except near the critical point where the free-volume
theory is less accurate. Lekkerkerker  _et al._ (1992); Dijkstra  _et al._
(2006) In Fig. 9.a, we show the $\eta_{\text{c}}$ profiles for $q=0.8$.
Surprisingly, we find weak oscillations at the "colloid-rich" liquid side of
the density profiles near the triple point. These oscillations were predicted
in fundamental measure theory, Brader  _et al._ (2001) but were not observed
in both a square-gradient density functional approach Brader and Evans (2000)
and in recent experiments on colloid-polymer mixtures. Vis  _et al._ (2020)
Here, we show for the first time such density oscillations, although weak, at
a free gas-liquid interface of a colloid-polymer mixture. The amplitude of
these oscillations are reduced significantly by thermally induced capillary-
wave fluctuations, Aarts, Schmidt, and Lekkerkerker (2004) making the
oscillations hard to detect in experiments Vis  _et al._ (2020) and
simulations. Chacón  _et al._ (2001) Such capillary wave fluctuations are
related to the ultra-low interfacial tensions of colloidal systems. Brader
_et al._ (2001)

To investigate if our ML potentials can also be used to measure the
interfacial gas-liquid tension $\gamma$, we use the so-called Test Area MC
technique (TAMC), which is based on thermodynamic perturbation theory. Gloor
_et al._ (2005); Vega and De Miguel (2007) In particular, we sample $\gamma$
in our MC simulations by performing test area perturbations once every cycle
and averaging over $5\times 10^{5}$ cycles. We report the interfacial tension
$\beta\gamma\sigma_{\text{c}}^{2}$ as a function of the difference of the
coexisting densities of the gas and liquid phases
$\eta_{\text{c}}^{\text{L}}-\eta_{\text{c}}^{\text{G}}$ in Fig. 9.b for a
colloid-polymer mixture with $q=0.8$. The reported error bars are obtained by
dividing the averaging run into 10 subaverages. We find qualitative agreement
with results for the true binary mixture reported by Vink et al. Verso  _et
al._ (2006)

##  V Conclusions

In conclusion, we have introduced a coarse-graining ML approach for colloidal
systems in which we trace out the degrees of freedom of the microscopic
species and fit the resulting effective many-body potential with a set of SFs
using simple linear regression. We have applied this approach to a model
suspension of colloidal particles and non-adsorbing polymer, and found that
the ML potential accurately describes the effective many-body potential for a
wide range of colloid densities and polymer fugacities. The ML potential can
therefore be assumed to be effectively state-independent and can even be used
to simulate direct phase coexistence. Given the computational efficiency of
the fitted model, we were able to test the validity of the ML approach by
measuring the gas-liquid binodals using MC direct-coexistence simulations. We
found good agreement with previous results obtained by using simulations of
the true binary mixture or using the full effective Hamiltonian that includes
all many-body interactions, even for the largest $q$ that we studied, where
four- and higher-body interactions are most pronounced. In addition, the
structure was also well-captured by the ML potentials, but deviations appear
when four- and higher-body contributions become predominant. Constructing
higher-body encoders (SFs) seems to be necessary to describe more accurately
the structure of the fluid. A generalization of the method to non-spherical
bodies is important to account for the effective many-body interactions of
anisotropic colloids, like mixtures of non-adsorbing polymers with colloidal
rods,Patti and Dijkstra (2009); Savenko and Dijkstra (2006) core-corona
nanorods,Campos-Villalobos, Dijkstra, and Patti (2021) polyhedral-shaped
particles,Henzie  _et al._ (2012) and superballs.Rossi  _et al._ (2015) It
will also be interesting to investigate whether this ML approach can be
extended to charged colloids, ligand-stabilized nanoparticles, starpolymers,
which will be subject of future work.

###### Acknowledgements.

G.C.V. acknowledges funding from The Netherlands Organisation for Scientific
Research (NWO) for the ENW PPS Fund 2018 – Technology Area Soft Advanced
Materials ENPPS.TA.018.002. M.D. has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (Grant agreement No. ERC-2019-ADG 884902 SoftML). L.F.
and E.B. acknowledge funding from The Netherlands Organisation for Scientific
Research (NWO) (Grant No. 16DDS004), and L.F. acknowledges funding from NWO
for a Vidi grant (Grant No. VI.VIDI.192.102).

## Conflict of interest

The authors have no conflicts to disclose.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

*

##  Appendix A Feature Selection Procedure

In this appendix we compare the efficiency of the feature selection procedure
used in this work (introduced in Ref. Boattini  _et al._ , 2020) with two of
the established selection methods by Imbalzano et al. Imbalzano  _et al._
(2018) In particular, we test the Pearson Correlation method (PC), and the
Farthest Point Sampling (FPS) scheme when using a linear regression scheme
like we did in this work.

The PC and FPS methods introduced by Imbalzano et al. are based solely on the
knowledge of the geometry of the particles’ environments, and do not rely on
the energy (or $V_{\text{f}}^{(3+)}$ in our case), nor on the performance of
the model that results from a given choice of the SFs. Instead, the common
idea behind those methods is to choose SFs which are as diverse as possible
by, e.g., minimizing their linear correlation or maximizing their difference
(in terms of their Euclidean distance), with the goal of minimizing the
redundancy in the selected subset. The method we use, instead, aims at
maximizing the correlation between the selected subset and the target
$V_{\text{f}}^{(3+)}$ by iteratively selecting SFs that add relevant
information (i.e. relevant to the task of fitting $V_{\text{f}}^{(3+)}$) to
the previously selected set, while penalizing both (i) highly correlated SFs
with only redundant information, and (ii) SFs which are sensitive to aspects
of the particle’ s environment that poorly correlate with
$V_{\text{f}}^{(3+)}$. The inclusion of point (ii) is arguably the main
difference with the other methods and guarantees a more efficient selection.

To demonstrate this, we show in Fig. 10, the performance of the three methods
in terms of the RMSE on the test set as a function of the number of selected
SFs for three different values of $q$. In all cases, the selection procedure
of Ref. Boattini  _et al._ , 2020 clearly outperforms the other methods,
leading to a higher accuracy of the fit. Note that, in all cases and for all
three methods, the first SF that is selected is the one with the largest
correlation with $V_{\text{f}}^{(3+)}$.

Figure 10: Root mean squared error RMSE as a function of the number of
selected SFs $N_{\text{s}}$ for a colloid-polymer mixture with size ratio
$q=0.4$, $0.8$ and $1.0$. We show the results for the selection procedure used
in our work (Ref. Boattini  _et al._ , 2020) and the PC and FPS schemes by
Imbalzano et al. Imbalzano  _et al._ (2018)

## References

  * Asakura and Oosawa (1954) S. Asakura and F. Oosawa, “On interaction between two bodies immersed in a solution of macromolecules,” J. Chem. Phys. 22, 1255–1256 (1954).
  * Gast, Hall, and Russel (1983) A. Gast, C. Hall, and W. Russel, “Polymer-induced phase separations in nonaqueous colloidal suspensions,” J. Colloid Interface Sci. 96, 251–267 (1983).
  * Meijer and Frenkel (1991) E. J. Meijer and D. Frenkel, “Computer simulation of polymer-induced clustering of colloids,” Phys. Rev. Lett. 67, 1110 (1991).
  * Lekkerkerker  _et al._ (1992) H. N. Lekkerkerker, W.-K. Poon, P. N. Pusey, A. Stroobants, and P. . Warren, “Phase behaviour of colloid+ polymer mixtures,” Europhys. Lett. 20, 559 (1992).
  * Ilett  _et al._ (1995) S. M. Ilett, A. Orrock, W. Poon, and P. Pusey, “Phase behavior of a model colloid-polymer mixture,” Phys. Rev. E 51, 1344 (1995).
  * Dijkstra, Brader, and Evans (1999) M. Dijkstra, J. M. Brader, and R. Evans, “Phase behaviour and structure of model colloid-polymer mixtures,” J. Phys. Condens. Matter 11, 10079 (1999).
  * Dijkstra (2001) M. Dijkstra, “Computer simulations of charge and steric stabilised colloidal suspensions,” Curr. Opin. Colloid Interface Sci. 6, 372–382 (2001).
  * Dijkstra and van Roij (1997) M. Dijkstra and R. van Roij, “Entropy-driven demixing in binary hard-core mixtures: From hard spherocylinders towards hard spheres,” Phys. Rev. E 56, 5594 (1997).
  * Buhot and Krauth (1998) A. Buhot and W. Krauth, “Numerical solution of hard-core mixtures,” Phys. Rev. Lett. 80, 3787 (1998).
  * Lobaskin and Linse (1999) V. Lobaskin and P. Linse, “Simulation of an asymmetric electrolyte with charge asymmetry 60: 1 using hard-sphere and soft-sphere models,” J. Chem. Phys. 111, 4300–4309 (1999).
  * Liu and Luijten (2004) J. Liu and E. Luijten, “Rejection-free geometric cluster algorithm for complex fluids,” Phys. Rev. Lett. 92, 035504 (2004).
  * Vink and Horbach (2004) R. Vink and J. Horbach, “Grand canonical monte carlo simulation of a model colloid–polymer mixture: Coexistence line, critical behavior, and interfacial tension,” J. Chem. Phys. 121, 3253–3258 (2004).
  * Bernard, Krauth, and Wilson (2009) E. P. Bernard, W. Krauth, and D. B. Wilson, “Event-chain monte carlo algorithms for hard-sphere systems,” Phys. Rev. E 80, 056704 (2009).
  * Dijkstra and Frenkel (1994) M. Dijkstra and D. Frenkel, “Evidence for entropy-driven demixing in hard-core fluids,” Phys. Rev. Lett. 72, 298 (1994).
  * Panagiotopoulos and Kumar (1999) A. Z. Panagiotopoulos and S. K. Kumar, “Large lattice discretization effects on the phase coexistence of ionic fluids,” Phys. Rev. Lett. 83, 2981 (1999).
  * Löwen, Madden, and Hansen (1992) H. Löwen, P. A. Madden, and J.-P. Hansen, “Ab initio description of counterion screening in colloidal suspensions,” Phys. Rev. Lett. 68, 1081 (1992).
  * Fushiki (1992) M. Fushiki, “Molecular-dynamics simulations for charged colloidal dispersions,” J. Chem. Phys. 97, 6700–6713 (1992).
  * Car and Parrinello (1985) R. Car and M. Parrinello, “Unified approach for molecular dynamics and density-functional theory,” Phys. Rev. Lett. 55, 2471 (1985).
  * Tata, Rajalakshmi, and Arora (1992) B. Tata, M. Rajalakshmi, and A. K. Arora, “Vapor-liquid condensation in charged colloidal suspensions,” Phys. Rev. Lett. 69, 3778 (1992).
  * Ito, Yoshida, and Ise (1994) K. Ito, H. Yoshida, and N. Ise, “Void structure in colloidal dispersions,” Science 263, 66–68 (1994).
  * Kepler and Fraden (1994) G. M. Kepler and S. Fraden, “Attractive potential between confined colloids at low ionic strength,” Phys. Rev. Lett. 73, 356 (1994).
  * Tata  _et al._ (1997) B. Tata, E. Yamahara, P. Rajamani, and N. Ise, “Amorphous clustering in highly charged dilute poly (chlorostyrene-styrene sulfonate) colloids,” Phys. Rev. Lett. 78, 2660 (1997).
  * Larsen and Grier (1997) A. E. Larsen and D. G. Grier, “Like-charge attractions in metastable colloidal crystallites,” Nature 385, 230–233 (1997).
  * Aarts, Schmidt, and Lekkerkerker (2004) D. G. Aarts, M. Schmidt, and H. N. Lekkerkerker, “Direct visual observation of thermal capillary waves,” Science 304, 847–850 (2004).
  * Vink, Horbach, and Binder (2005) R. Vink, J. Horbach, and K. Binder, “Capillary waves in a colloid-polymer interface,” J. Chem. Phys. 122, 134905 (2005).
  * Moussaid  _et al._ (1999) A. Moussaid, W. Poon, P. Pusey, and M. Soliva, “Structure of marginal and fully developed colloidal liquids,” Phys. Rev. Lett. 82, 225 (1999).
  * Vis  _et al._ (2020) M. Vis, K. J. Brouwer, A. Gonzalez Garcia, A. V. Petukhov, O. Konovalov, and R. Tuinier, “Quantification of the structure of colloidal gas–liquid interfaces,” J. Phys. Chem. Lett. 11, 8372–8377 (2020).
  * Dijkstra  _et al._ (2006) M. Dijkstra, R. van Roij, R. Roth, and A. Fortini, “Effect of many-body interactions on the bulk and interfacial phase behavior of a model colloid-polymer mixture,” Phys. Rev. E 73, 041404 (2006).
  * Kobayashi  _et al._ (2019) H. Kobayashi, P. B. Rohrbach, R. Scheichl, N. B. Wilding, and R. L. Jack, “Correction of coarse-graining errors by a two-level method: Application to the asakura-oosawa model,” J. Chem. Phys. 151, 144108 (2019).
  * Rupp  _et al._ (2012) M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. Von Lilienfeld, “Fast and accurate modeling of molecular atomization energies with machine learning,” Phys. Rev. Lett. 108, 058301 (2012).
  * Behler (2016) J. Behler, “Perspective: Machine learning potentials for atomistic simulations,” J. Chem. Phys. 145, 170901 (2016).
  * Musil  _et al._ (2021) F. Musil, A. Grisafi, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, “Physics-inspired structural representations for molecules and materials,” arXiv preprint arXiv:2101.04673 (2021).
  * Verso  _et al._ (2006) F. L. Verso, R. Vink, D. Pini, and L. Reatto, “Critical behavior in colloid-polymer mixtures: Theory and simulation,” Phys. Rev. E 73, 061407 (2006).
  * Behler and Parrinello (2007) J. Behler and M. Parrinello, “Generalized neural-network representation of high-dimensional potential-energy surfaces,” Phys. Rev. Lett. 98, 146401 (2007).
  * Boattini  _et al._ (2020) E. Boattini, N. Bezem, S. N. Punnathanam, F. Smallenburg, and L. Filion, “Modeling of many-body interactions between elastic spheres through symmetry functions,” J. Chem. Phys. 153, 064902 (2020).
  * Asakura and Oosawa (1958) S. Asakura and F. Oosawa, “Interaction between particles suspended in solutions of macromolecules,” J. Polym. Sci. 33, 183–192 (1958).
  * Vrij (1976) A. Vrij, “Polymers at interfaces and the interactions in colloidal dispersions,” Pure Appl. Chem. 48, 471–483 (1976).
  * Dijkstra and van Roij (2002) M. Dijkstra and R. van Roij, “Entropic wetting and many-body induced layering in a model colloid-polymer mixture,” Phys. Rev. Lett. 89, 208303 (2002).
  * Behler (2011) J. Behler, “Atom-centered symmetry functions for constructing high-dimensional neural network potentials,” J. Chem. Phys. 134, 074106 (2011).
  * Note (1) The configurations we use are spaced far enough apart in number of MC cycles that the particles have had a chance to diffuse significantly between successive snapshots. Even in the high density fluid phase at $\eta_{\text{c}}=0.40$, the particles travel more than 6 times their diameters between successive configurations.
  * Behler (2015) J. Behler, “Constructing high-dimensional neural network potentials: A tutorial review,” Int. J. Quantum Chem. 115, 1032–1050 (2015).
  * Khaliullin  _et al._ (2010) R. Z. Khaliullin, H. Eshet, T. D. Kühne, J. Behler, and M. Parrinello, “Graphite-diamond phase coexistence study employing a neural-network mapping of the ab initio potential energy surface,” Phys. Rev. B 81, 100103 (2010).
  * Eshet  _et al._ (2012) H. Eshet, R. Z. Khaliullin, T. D. Kühne, J. Behler, and M. Parrinello, “Microscopic origins of the anomalous melting behavior of sodium under high pressure,” Phys. Rev. Lett. 108, 115701 (2012).
  * Kapil, Behler, and Ceriotti (2016) V. Kapil, J. Behler, and M. Ceriotti, “High order path integrals made easy,” J. Chem. Phys. 145, 234103 (2016).
  * Cheng, Behler, and Ceriotti (2016) B. Cheng, J. Behler, and M. Ceriotti, “Nuclear quantum effects in water at the triple point: Using theory as a link between experiments,” J. Phys. Chem. Lett. 7, 2210–2215 (2016).
  * Imbalzano  _et al._ (2018) G. Imbalzano, A. Anelli, D. Giofré, S. Klees, J. Behler, and M. Ceriotti, “Automatic selection of atomic fingerprints and reference configurations for machine-learning potentials,” J. Chem. Phys. 148, 241730 (2018).
  * Dijkstra, van Roij, and Evans (2000) M. Dijkstra, R. van Roij, and R. Evans, “Effective interactions, structure, and isothermal compressibility of colloidal suspensions,” J. Chem. Phys. 113, 4799–4807 (2000).
  * Louis (2002) A. Louis, “Beware of density dependent pair potentials,” J. Phys. Condens. Matter 14, 9187 (2002).
  * Brader  _et al._ (2001) J. Brader, R. Evans, M. Schmidt, and H. Löwen, “Entropic wetting and the fluid-fluid interface of a model colloid-polymer mixture,” J. Phys. Condens. Matter 14, L1 (2001).
  * Brader and Evans (2000) J. Brader and R. Evans, “The fluid-fluid interface of a model colloid-polymer mixture,” EPL (Europhysics Letters) 49, 678 (2000).
  * Chacón  _et al._ (2001) E. Chacón, M. Reinaldo-Falagán, E. Velasco, and P. Tarazona, “Layering at free liquid surfaces,” Phys. Rev. Lett. 87, 166101 (2001).
  * Gloor  _et al._ (2005) G. J. Gloor, G. Jackson, F. J. Blas, and E. de Miguel, “Test-area simulation method for the direct determination of the interfacial tension of systems with continuous or discontinuous potentials,” J. Chem. Phys. 123, 134703 (2005).
  * Vega and De Miguel (2007) C. Vega and E. De Miguel, “Surface tension of the most popular models of water by using the test-area simulation method,” J. Chem. Phys. 126, 154707 (2007).
  * Patti and Dijkstra (2009) A. Patti and M. Dijkstra, “Do multilayer crystals nucleate in suspensions of colloidal rods?” Phys. Rev. Lett. 102, 128301 (2009).
  * Savenko and Dijkstra (2006) S. Savenko and M. Dijkstra, “Phase behavior of a suspension of colloidal hard rods and nonadsorbing polymer,” J. Chem. Phys. 124, 234902 (2006).
  * Campos-Villalobos, Dijkstra, and Patti (2021) G. Campos-Villalobos, M. Dijkstra, and A. Patti, “Nonconventional phases of colloidal nanorods with a soft corona,” Phys. Rev. Lett. 126, 158001 (2021).
  * Henzie  _et al._ (2012) J. Henzie, M. Grünwald, A. Widmer-Cooper, P. L. Geissler, and P. Yang, “Self-assembly of uniform polyhedral silver nanocrystals into densest packings and exotic superlattices,” Nat. Mater. 11, 131–137 (2012).
  * Rossi  _et al._ (2015) L. Rossi, V. Soni, D. J. Ashton, D. J. Pine, A. P. Philipse, P. M. Chaikin, M. Dijkstra, S. Sacanna, and W. T. Irvine, “Shape-sensitive crystallization in colloidal superball fluids,” Proc. Natl. Acad. Sci. U.S.A. 112, 5286–5290 (2015).