THE JOURNAL OF CHEMICAL PHYSICS 143, 084501 (2015)
Pattern formation in binary ﬂuid mixtures induced by shortrange
competing interactions
Cecilia Bores,
1
Enrique Lomba,
1,2,a)
Aurélien Perera,
2
and Noé G. Almarza
1
1
Instituto de Química Física Rocasolano, CSIC, Serrano 119, E28006 Madrid, Spain
2
Laboratoire de Physique Théorique de la Matiére Condensée (UMR CNRS 7600),
Université Pierre et Marie Curie, 4 Place Jussieu, F75252 Paris Cedex 05, France
(Received 29 June 2015; accepted 23 July 2015; published online 24 August 2015)
Molecular dynamics simulations and integral equation calculations of a simple equimolar mixture of
diatomic molecules and monomers interacting via attractive and repulsive shortrange potentials show
the existence of pattern formation (microheterogeneity), mostly due to depletion forces away from
the demixing region. Eﬀective sitesite potentials extracted from the pair correlation functions using
an inverse Monte Carlo approach and an integral equation inversion procedure exhibit the features
characteristic of a shortrange attractive and a longrange repulsive potential. When charges are
incorporated into the model, this becomes a coarse grained representation of a room temperature ionic
liquid, and as expected, intermediate range order becomes more pronounced and stable.
C
2015 AIP
Publishing LLC. [http://dx.doi.org/10.1063/1.4928524]
I. INTRODUCTION
Spontaneous pattern formation is a feature present in a
diverse collection of physical, chemical, and biological sys
tems.
1
In spite of the diverse nature of these systems, the
appearance of the emerging microphases is quite similar: in
2D systems, circular droplets, stripes, or “bubbles” occur, and
in 3D systems, one may ﬁnd spherical droplets, sheets, or
tubes. In some cases, the patterns appear as transient states
due to energy or mass ﬂuctuations that occur in the process
of spinodal decomposition, but sometimes, these states can be
stabilized due to the presence of competitive interactions, in
which one of the interactions is responsible for inhibiting the
phase separation.
2,3
The understanding of this selforganizing capability of
soft and ﬂuid matter is critical for a wide panoply of appli
cations of great relevance nowadays. These selfassembly
mechanisms play a crucial role in processes involving protein
solutions in food products,
4,5
therapeutic monoclonal anti
bodies,
6–8
nanolithography,
9
or gelation processes.
10
In the realm of colloidal science, systems with extremely
short ranged repulsive interactions are often used as an exper
imental realization of the hard sphere ﬂuid,
11
a system noto
rious for its theoretical interest. On the other hand, the addi
tion of nonadsorbing polymers to the colloidal solution typi
cally activates an attractive interparticle interaction, due to the
depletion mechanism. Moreover, changing the concentration
and molecular weight of the polymer, the attraction range
and strength of the colloidcolloid interaction can be tuned.
Clustering is to be expected due to the presence of the attractive
forces,
12,13
but in principle, it would correspond to metastable
states and/or irreversible processes of kinetic nature. Never
theless, microphases formed by clusters and percolating struc
tures can be stabilized in protein solutions and colloidpolymer
a)
enrique.lomba@csic.es
mixtures both in experiment
14
and in theoretical descriptions
15
due to the presence of additional repulsive interactions stem
ming from electrostatic forces. An extreme example of the
stabilizing role of charges is the nanostructural organization
that appears in room temperature ionic liquids (RTILs).
16
In
fact, it has been shown that long range repulsive interactions
alone can give rise to nanostructural order,
17
the driving force
of attractive interactions to induce spontaneous aggregation
being replaced by external forces (e.g., pressure).
In the caseofcolloidal systems, in which chargedcolloidal
particles are screened by ions in the solvent, the colloid
colloid interaction has been shown on theoretical grounds to be
adequately represented by a Yukawa potential
18,19
according
to the DerjaguinLandauVerweyOverbeek (DLVO) theory.
Following this, numerous works have resorted to potentials
with a combination of a short range attraction and a long range
repulsion (SALR) in the form of a double Yukawa,
20,21
or a
LennardJones (LJ) plus a Yukawa interaction
2,22,23
in order to
model the spontaneous emergence of microstructured patterns
in ﬂuids. On the other hand, back in 1999, Sear et al.
24
made
use of an empirical two exponential form with SALR char
acteristics in order to explain the experimental appearance of
stable microphases of nanoparticles at the airwater interface.
This potential has been studied in depth in model systems, both
in bulk and in conﬁnement,
25–29
and as a rough approximation
to account for vegetation patterns in ecosystems with limited
resources.
30
In this work, we will explore the possibility of pattern
formation in a system in which only short ranged forces are
present. Our model system, composed of heteronuclear dimers
and monomers, combines attractive and repulsive potentials
so as to mimic the interactions present in RTILs, but without
electrostatic forces. To that aim, we have performed extensive
molecular dynamics simulations in the canonical and in the
isothermalisobaric ensembles. We will address the emergence
of intermediate range order (IRO)analyzing the behaviorof the
00219606/2015/143(8)/084501/11/$30.00 143, 0845011 © 2015 AIP Publishing LLC
Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rightsandpermissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49
0845012 Bores et al. J. Chem. Phys. 143, 084501 (2015)
partial and concentrationconcentration structure factors and
performing a cluster analysis for various degrees of asymmetry
in the sites of the diatomic particles. Reference Interaction
Site Model (RISM) integral equation calculations have also
been carried out and are shown to agree remarkably well with
the simulation results. By means of an Inverse Monte Carlo
(IMC) approach,
31
we have extracted eﬀective interactions
from the pair correlation functions of the simulated mixtures.
For comparison, another set of eﬀective potentials has been
obtained from the RISM results using an integral equation
inversion procedure. We will see that despite the fact that all
interactions at play are short ranged, their net eﬀect leading to
the pattern formation (microheterogeneity, or microstructure
segregation at the nanoscale) translates into the appearance of
eﬀective interactions that agree with the characteristic trends
of a short range attraction and a long range repulsion, i.e., a
SALR potential. We have found that the eﬀective potentials ex
tracted from the simulation and those derived by the theoretical
approach agree remarkably well. Finally, we have analyzed the
role of charges on our model, which in fact by the addition of
electrostatic sitesite interactions becomes a rough representa
tion of a RTIL. As expected, charges will be shown to enhance
the pattern formation and the stability of the nanostructured
phases.
The rest of the paper can be sketched as follows. In Sec. II
we introduce the model in full detail and brieﬂy summarize the
methodology. In Section III, we introduce our most signiﬁcant
results. Conclusions and future prospects are to be found in
Section IV.
II. MODEL AND METHODS
Our model consists in an equimolar ﬂuid mixture of two
diﬀerent species, a twosite dimer AB and a monomer C. The
dimers are represented by a two center LJ sitesite potential, in
which the sites are separated by a distance l. Our monomers
also interact via LJ potentials. In all cases, the interactions are
cut and shifted at a distance r
c
, by which the explicit form of
the sitesite potentials is
u
i j
(r) = 4ϵ
σ
i j
r
12
−
σ
i j
r
6
−
σ
i j
r
c
12
+
σ
i j
r
c
6
if r < r
c
, (1)
and u
i j
(r) = 0 otherwise. Our model is to a certain degree
inspired by the simple coarsegrained model for imidazolium
based RTIL of Merlet et al.
32
We will see to what extent a
simple model, with just two sites and purely short ranged inter
actions can reproduce the presence of nanostructural order as
found in RTILs. To that aim, we will however preserve the
attractive/repulsive character of the interactions in the RTIL. In
our model then, C monomers would correspond to anions, AB
dimers to the molecular cations, with the imidazolium ring that
contains the positive charge, being represented by site A, and
the nonpolar tail, by the larger site B. This implies that AA and
CC interactions will be repulsive, BB and AC are attractive,
ﬁnally BC and AB interactions are also repulsive. For the
sizes of A and C particles, we have chosen σ
AA
= σ
CC
= 4 Å,
the elongation of the dimer l = 8 Å. The AB distances of the
TABLE I. LennardJones potential parameters.
Particle
i
Particle
j Interaction ϵ (kJ/mol) σ
i j
r
c
A A Repulsive 2.092 4.0 Å 2
1/6
· σ
A A
A B Repulsive 2.092 (σ
A A
+ σ
B B
)/2 2
1/6
· σ
AB
A C Attractive 2.092 4.0 Å 3 · σ
B B
B B Attractive 2.092 σ
B B
3 · σ
B B
B C Repulsive 2.092 (σ
B B
+ σ
CC
)/2 2
1/6
· σ
BC
C C Repulsive 2.092 4.0 Å 2
1/6
· σ
CC
dimers are ﬁxed as constraints of the equations of motion. The
LJ well is set to ϵ = 2.092 kJ/mol, identical for all interactions.
Since the size of the nonpolar tail is essential to determine the
nanostructural ordering,
16
we have considered various sizes
for σ
BB
(with σ
BB
> σ
AA
always). For the attractive inter
actions, we have truncated and shifted the LJ potential at r
c
= 3σ
BB
. For the repulsive interactions, we have simply used r
c
= 2
1/6
σ
i j
, thus deﬁning purely repulsive soft spheres follow
ing the prescription of Weeks, Chandler, and Andersen
(WCA).
33
The complete set of parameters for all interactions
is summarized in Table I. Finally, in order to analyze the eﬀect
of charges on the intermediate range order, we have considered
explicitly the same model with a positive charge +q on the A
sites and a corresponding negative charge −q on the monomers.
The value of q is varied between 0 and 0.25e, where e is the
elementary electron charge. Again these values are of the same
order as those considered in the model of Ref. 32.
A. Simulations and analysis
We have carried out extensive molecular dynamics simu
lations of the system previously described using the LAMMPS
package,
34–36
in the canonical and isothermalisobaric ensem
bles using a NoseHoover thermostat and barostat.
37
Our sam
ples contained 16 384 particles (samples of up to 65 536 parti
cles were investigated and no signiﬁcant size dependence was
found). For simplicity, we considered equal masses for the
three interaction centers: m
A
= m
B
= m
c
= 16 g mol
−1
. Initial
thermalization runs at a temperature of 226 K were 2 × 10
6
steps long, with a time step of 1 fs. Production runs were
5 × 10
6
steps long, and averages were carried out every 5000
steps.
One of the problems one can encounter when performing
canonical simulations in this type of system is the occur
rence of phase transitions, either vaporliquid equilibria or
demixing. In order to guarantee that the states under consid
eration correspond to thermodynamic equilibrium conditions,
and consequently, any potential intermediate range order is not
the result of a spinodal decomposition, we have run additional
isothermalisobaric simulations and analyzed the volume ﬂuc
tuation of the samples. In this way, one can avoid those states
that lie inside the liquidvapor spinodal. Moreover, one can
compute the partial structure factors, deﬁned as
S
i j
(k) = x
i
δ
i j
+ x
i
x
j
ρ
g
i j
(r) − 1
e
−kr
dr, (2)
where ρ is the total number density, δ
i j
is a Kronecker δ, and
x
i
is the molar fraction of component i. Here, sites A and B
are considered as diﬀerent particles and g
i j
is the atomatom
Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rightsandpermissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49
0845013 Bores et al. J. Chem. Phys. 143, 084501 (2015)
pair distribution function. Our samples are large enough to
allow for an accurate integration of the pair distribution func
tions, and the results are consistent with direct ksampling.
Notice that as far as Eq. (2) is concerned, x
A
= x
B
= x
C
= 1/3;
hence, in the large k limit, all structure factors will tend to
1/3. From the partial structure factors, it is possible to evaluate
the concentrationconcentration structure factor introduced by
Bathia and Thornton,
38
for which we have deﬁned
S
cc
(k) = x
2
AB
S
CC
(k) + x
2
C
S
AB− AB
(k) − 2x
AB
x
C
S
C− AB
(k),
(3)
where now one has to consider explicitly the structure fac
tors corresponding to the molecular species AB, and as a
consequence, x
C
= x
AB
= 1/2. We can simply approximate
g
AB− AB
= g
BB
and g
C− AB
= g
C B
, as if the scattering length
or form factor of A sites was negligible compared to that of
B sites. This is in principle not unreasonable given the much
larger size of the B sites, but in a realistic situation, one should
take explicitly into account the true scattering lengths or form
factors of sites A and B. Now, one has to correct for the diﬀerent
values of the molar fraction when AB is considered as a single
species and Eq. (2) is used in (3). In this way, lim
k→ ∞
S
cc
(k)
= x
c
x
AB
= 1/4. With all this in mind, the presence of a diver
gence when k → 0 in S
cc
(k) is a signal of a demixing transi
tion, so this quantity will be essential to assess the stability of
the thermodynamic states chosen for our simulations.
Finally, back to the vaporliquid transition, one can ana
lyze the corresponding kdependent linear response suscepti
bility in density ﬂuctuations, namely,
39
ρk
B
T χ
T
(k) =
S(k)
ab
(x
a
x
b
)S(k)
ab
, (4)
whose k = 0 limit is precisely the isothermal compressibility.
In Eq. (4), k
B
is Boltzmann’s constant, T the absolute temper
ature, and the elements of the matrix S
i j
are just the partial
structure factors as deﬁned in Eq. (2).  . . .  denotes the matrix
determinant and  . . . 
ab
the corresponding minor of the matrix
S(k). The presence of a divergence—or a substantial increase
in χ
T
(0)—is a clear indication of the vicinity of a vaporliquid
transition. A careful monitoring of this quantity together with
the use of NPT simulations provides a reliable assessment of
the stability of the state points under consideration during the
simulation runs.
All systems and conditions studied in this work are sum
marized in Table II. In the case of system 8, when increasing
the charge from 0.10e to 0.25e, the conditions of temperature
and density corresponding to systems 3, 6, and 7 lie in the
twophase region. Consequently, we resorted to an isothermal
isobaric simulation at low positive pressure to achieve thermo
dynamic equilibrium conditions in our system with q = 0.25e.
The ﬁnal value of the total particle density achieved in this way
is indicated in Table II.
B. Inverse Monte Carlo method
With the pair correlation functions produced along the
simulation runs and the corresponding statistical uncertainties
calculated using block averages, we have used the IMC proce
dure proposed by Almarza and Lomba
31
in order to produce
TABLE II. Potential parameters and thermodynamic state variables for the
systems under study.
Potential thermodynamic state
q(e) σ
B
(Å) ρ (Å
−3
) T (K) P (MPa)
System 1 0 8.0 0.001 226.4 27.05
System 2 0 8.0 0.001 25 226.5 39.5
System 3 0 8.0 0.001 5 226.5 59.4
System 4 0 9.0 0.001 226.5 30.4
System 5 0 9.0 0.001 25 226.4 53.2
System 6 0 9.0 0.001 5 226.4 96.7
System 7 0.1 8.0 0.001 5 226.4 39.4
System 8 0.25 8.0 0.001 95 226.3 0.61
single component sitesite eﬀective potentials able to repro
duce the microscopic structure exhibited by our mixture model.
The procedure starts from a simple approximation βu
eﬀ
in
(r)
= − log g(r) and proceeds to modify the pair potential along the
simulation run in such a way that the calculated g
eﬀ
(r) matches
the input g(r). Explicit details of the method can be found in
Ref. 31. In our case, we have used a total of 4000 particles.
The procedure of inversion was carried out in 20 stages. In
the last stages, the eﬀective potentials hardly varied, and the
convergence between input and calculated g(r)’s according to
the prescription of Ref. 31 was achieved successfully in all the
cases.
In this way, one can use as input of the IMC procedure
either g
AA
(r), g
BB
(r), or g
CC
(r) and obtain a corresponding
set of u
eﬀ
AA
(r), u
eﬀ
BB
(r), and u
eﬀ
CC
(r), which will obviously be
diﬀerent, but in the case of emergence of intermediate range,
order should exhibit some common features.
C. RISM integral equation
The sitesite correlations are obtained by solving the usual
set of 2 equations, the sitesite OrnsteinZernike (SSOZ) equa
tion and the closure equation, which we choose here to be the
sitesite hypernetted (SSHNC) equation. The SSOZ equation
for the present system is explicitly given in the matrix form
(W +
ρ
3
H)(W
−1
−
ρ
3
C) = I, (5)
where the 3 × 3 matrix H (or C) has for elements H
i j
=
˜
h
i j
(k)
(or C
i j
= ˜c
i j
(k)), the Fourier transform (FT) of the sitesite pair
correlation functions h
i j
(r) = g
i j
(r) − 1 (or the direct correla
tion function c
i j
(r)), where the indices i and j stand for one
of the sites A, B, and C. The matrix W represents the intra
molecular correlations, which for the present system gives
W =
*
.
.
.
,
˜w
AA
˜w
AB
˜w
AC
˜w
AB
˜w
BB
˜w
BC
˜w
AC
˜w
BC
˜w
CC
+
/
/
/

=
*
.
.
.
,
1 j
0
(kl) 0
j
0
(kl) 1 0
0 0 1
+
/
/
/

, (6)
where j
0
(x) is a spherical Bessel function. The matrix I is the
identity matrix. The SSHNC equations are written as
g
i j
(r) = exp
−
u
i j
(r)
k
B
T
+ h
i j
(r) − c
i j
(r)
, (7)
and there are 9 such independent equations to solve.
Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rightsandpermissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49
0845014 Bores et al. J. Chem. Phys. 143, 084501 (2015)
Both equations are approximate and their respective
inconsistencies have been discussed many times in the past
literature.
39,40
Based on empirical evidence from the literature,
we expect that the correlations obtained through these equa
tions for the present systems, both charged and uncharged,
should be relatively good for the short range part, but perhaps
not at long range. We are particularly interested to see if the
correlations related to the appearance of the local structures
can be reproduced by this theory. The structure factor deﬁned
in Eq. (3) is the appropriate function for this purpose, as illus
trated in Sec. III.
The practical solution of these equations consists in dis
cretizing all the functions on an equidistant grid, both in r and
k space. We use 2048 points with a rgrid of ∆r = 0.01σ
A
,
which is enough for the present case to properly describe the
asymptotic behavior of the correlations in direct and reciprocal
space. The set of two equations are solved iteratively following
techniques well documented in the literature.
It is also possible to obtain the eﬀective potentials which
would correspond to the equivalent onecomponent represen
tation of the system. This is achieved by imposing the pair
correlation function to be the desired sitesite correlation,
namely, g(r) = g
X X
(r), in the set of the two integral equations
for the 1component system and solve these equations for the
direct correlation function and eﬀective pair interaction. The
direct correlation function can be obtained through the OZ
equation for 1component system (which is an exact relation),
(1 + ρ
S
˜
h(k))(1 − ρ
S
˜c(k)) = 1, (8)
where h(r) = g(r) − 1 = h
X X
(r) = g
X X
(r) − 1, and the den
sity ρ
S
is that of the eﬀective 1 component made solely of
sites X . Once c(r) is obtained, one solves the HNC closure,
which has the same form as Eq. (7), but now for the eﬀective
interaction u
e f f
(r), one gets
u
e f f
(r) = −k
B
T
[
ln g
X X
(r) + h
X X
(r) − c(r)
]
. (9)
III. RESULTS
A. Pair structure
Here, we have analyzed the eﬀect of the molecular geom
etry on the nanostructure formation changing the diameter of
σ
BB
. We have ﬁrst considered σ
BB
= 8 Å, 9 Å, 10 Å, and 12 Å.
Some snapshots of conﬁgurations for varying σ
BB
are depicted
in Figure 1. We have found that for σ
BB
> 9 Å, clustering
or microheterogeneity of C particles can only be appreciated
when the packing of the B sites is so high that it resembles
that of a solid. In fact in this case, the height of the ﬁrst peak
of S
BB
(k) exceeds 2.7, which according to the HansenVerlet
rule
41
indicates that freezing conditions have been reached.
Moreover, the prepeak in the structure factor characteristic of
the presence of IRO is absent from S
BB
(k). The clustering of
C particles results from a merely steric eﬀect, since these are
restricted to occupy the holes between the large B particles.
These eﬀects can be appreciated in the snapshots of Figure 1,
where the dense packing of B sites (red spheres) is readily
apparent.
For the reason mentioned above, we will concentrate on
the results for σ
BB
= 8 Å and 9 Å. Already in the correspond
ing snapshot of Figure 1, one can appreciate the formation of a
bicontinuous network of percolating clusters, connecting both
AB dimers and C monomers. By bicontinuous network, we
mean that the clusters formed by Bsites and C particles will
be seen to both span practically the whole sample, forming two
continuous interpenetrated percolating microphases. This can
be analyzed from a more quantitative perspective by ﬁrst taking
a look at the corresponding pair distribution functions and par
tial structure factors, which are depicted in Figures 2 and 3,
respectively, for Systems 1–6. Focusing ﬁrst on the g
AA
pair
distribution function, one ﬁrst appreciates the large exclusion
hole after the ﬁrst layer, which is a simple consequence of the
large size of Bsites. Obviously, the exclusion hole grows with
the size of the Bsites, as can be seen when comparing ﬁgures
FIG. 1. Snapshots of conﬁgurations for total particle density ρ = 0.001 25 Å
−3
and temperature T = 226.45 K for two Bsite diameters. As the size of
Bsites grows, C monomers cluster in the cavities formed by the Bsites due to excluded volume eﬀects. All other diameters and total density are kept ﬁxed,
σ
A A
= σ
CC
= 4.0 Å. (a) σ
B B
= 8 Å. (b) σ
B B
= 12 Å.
Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rightsandpermissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49
0845015 Bores et al. J. Chem. Phys. 143, 084501 (2015)
FIG. 2. The ﬁgures show the radial
distribution functions for A, B, and C
particles, respectively. Column (a) cor
responds to σ
B B
= 8 Å for system
3 (theory vs. simulation) and column
(b) presents the simulations results for
systems 4–6 for σ
B B
= 9 Å. Total den
sity is indicated in the legend. Simu
lation results are represented by solid
lines and dasheddotted curves corre
spond to integral equation calculations.
on the left and right columns. Correlations between Asites
extend up to ﬁve σ
AA
, and the width of the g
CC
correlation
is ≈2σ
CC
. These features hint at the presence of some degree
of IRO. BB correlations (graphs in the middle row) behave
like those of a dense ﬂuid, and no apparent sign of clustering
or IRO is evident. In contrast, the wide ﬁrst peak of g
CC
is
characteristic of clusters of particles conﬁned in cavities, in
this case formed by Bsites. This eﬀect, as mentioned before,
is maximized for the largest σ
BB
. We will see later that these
clusters of partly occluded Cparticles are connected, forming
a three dimensional percolating structure.
If we take now a look at the partial structure factors, we
immediately appreciate a feature characteristic of the emer
gence of IRO, namely, the presence of a prepeak at 0.25 Å
−1
.
This corresponds to correlations in the range of 25 Å, the
distance at which any sign of structure of the pair distribution
function dies out. Interestingly, the prepeak is almost absent in
S
AA
, except for a small maximum visible for the σ
BB
= 9 Å
and the highest density. This quantity shows otherwise very
little structure for k > 0.5 Å
−1
. As seen in the g
AA
’s, the
most relevant feature in the AA correlations is the exclu
sion hole due to the presence of the Bsites. In contrast, S
BB
does exhibit a prepeak, even when no evidence of IRO was
visible in g
BB
. This prepeak is more apparent in the mono
mer structure factor S
CC
. When the density is lowered, the
prepeak in the Bsite structure factor shifts to lower kvalues
and vanishes at ρ = 0.001 Å
−3
. In the case of S
CC
, the posi
tion of the prepeak is preserved, but its magnitude decreases.
In Figure 4, the corresponding concentrationconcentration
structure factor is displayed. The prepeak at k
0
≈ 0.25 Å
−1
is
preserved, although its magnitude decreases when the total
density is lowered. In contrast, no increase when k → 0 is
visible. This implies that we are encountering concentration
ﬂuctuations inducing spatial inhomogeneities, but no demix
ing transition. In Figure 5, we have plotted the kdependent
susceptibility corresponding to density ﬂuctuations. The pre
peak is visible except for the lowest density, which implies
that density inhomogeneities with a spatial patterns are also
correlated with the corresponding concentration inhomogene
ities. But now, the k → 0 behavior is diﬀerent. As density is
decreased, the susceptibility (i.e., the isothermal compress
ibility) grows, an indication of the vicinity of a vaporliquid
transition. This means, that lowering the density from the value
of ρ = 0.001 Å
−1
at the same temperature could move the sys
tem across the spinodal curve into the twophase region. Our
analysis indicates that the thermodynamic conditions we have
simulated can be considered equilibrium states. Moreover, we
have conﬁrmed that the results do not have a signiﬁcant sample
size dependence, by which metastability can also be ruled out.
The sitesite correlation functions and structure factors
obtained from the RISM theory are represented in dashed
lines in Figs. 23. It is seen that the agreement is excellent
in most cases, particularly in what concerns the BB and CC
correlations. The AA correlations are systematically underes
timated near contact and overestimated at larger distances. The
most signiﬁcant diﬀerences are seen for the structure factors
in Fig. 3. Integral equations tend to exaggerate concentration
ﬂuctuations and often tend to interpret small aggregate forma
tions as such.
42,43
We observe here a similar trend for the low
density case ρ = 0.001 Å
−3
, for which ﬂuctuations compete the
most with aggregate formation. The prediction of aggregation,
through the prepeak is in very good agreement with simula
tions for the highest density ρ = 0.0015 Å
−3
, precisely when
the denser packing tends to favor aggregation. This is also in
line with previous observations of similar type of behavior for
model ionic liquids. These features are a direct consequence of
the fact that the HNC closure approximation misses high order
correlations, hence high order cluster contributions, which are
represented in the bridge term b
i j
(r) that is neglected in the
exponential of Eq. (7). We observe that in all cases, the k = 0
Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rightsandpermissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49