scispace - formally typeset
Search or ask a question
Journal ArticleDOI

Pattern formation in binary fluid mixtures induced by short-range competing interactions

24 Aug 2015-Journal of Chemical Physics (AIP Publishing)-Vol. 143, Iss: 8, pp 084501-084501
TL;DR: Molecular dynamics simulations and integral equation calculations of a simple equimolar mixture of diatomic molecules and monomers interacting via attractive and repulsive short-range potentials show the existence of pattern formation (microheterogeneity), mostly due to depletion forces away from the demixing region.
Abstract: Molecular dynamics simulations and integral equation calculations of a simple equimolar mixture of diatomic molecules and monomers interacting via attractive and repulsive short-range potentials show the existence of pattern formation (microheterogeneity), mostly due to depletion forces away from the demixing region. Effective site-site 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 short-range attractive and a long-range 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.

Summary (3 min read)

Introduction

  • Effective site-site 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 short-range attractive and a long-range 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.
  • The authors have found that the effective potentials extracted from the simulation and those derived by the theoretical approach agree remarkably well.

II. MODEL AND METHODS

  • The authors monomers also interact via LJ potentials.
  • The authors will see to what extent a simple model, with just two sites and purely short ranged interactions can reproduce the presence of nano-structural order as found in RTILs.
  • To that aim, the authors will however preserve the attractive/repulsive character of the interactions in the RTIL.
  • In their 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 non-polar tail, by the larger site B.
  • Again these values are of the same order as those considered in the model of Ref. 32.

A. Simulations and analysis

  • The authors have carried out extensive molecular dynamics simulations of the system previously described using the LAMMPS package,34–36 in the canonical and isothermal-isobaric ensembles using a Nose-Hoover thermostat and barostat.37.
  • One of the problems one can encounter when performing canonical simulations in this type of system is the occurrence of phase transitions, either vapor-liquid equilibria or demixing.
  • The final value of the total particle density achieved in this way is indicated in Table II.
  • With the pair correlation functions produced along the simulation runs and the corresponding statistical uncertainties calculated using block averages, the authors have used the IMC procedure proposed by Almarza and Lomba31 in order to produce single component site-site effective potentials able to reproduce the microscopic structure exhibited by their mixture model.

C. RISM integral equation

  • The site-site correlations are obtained by solving the usual set of 2 equations, the site-site Ornstein-Zernike (SSOZ) equation and the closure equation, which the authors choose here to be the site-site hypernetted (SS-HNC) equation.
  • Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions.
  • Based on empirical evidence from the literature, the authors expect that the correlations obtained through these equations for the present systems, both charged and uncharged, should be relatively good for the short range part, but perhaps not at long range.
  • This is achieved by imposing the pair correlation function to be the desired site-site correlation, namely, g(r) = gXX(r), in the set of the two integral equations for the 1-component system and solve these equations for the direct correlation function and effective pair interaction.

A. Pair structure

  • Here, the authors have analyzed the effect of the molecular geometry on the nanostructure formation changing the diameter of σBB.
  • Interestingly, the prepeak is almost absent in SAA, except for a small maximum visible for the σBB = 9 Å and the highest density.
  • As seen in the gAA’s, the most relevant feature in the AA correlations is the exclusion hole due to the presence of the B-sites.
  • When the density is lowered, the prepeak in the B-site structure factor shifts to lower k-values and vanishes at ρ = 0.001 Å−3.
  • The authors observe here a similar trend for the low density case ρ = 0.001 Å−3, for which fluctuations compete the most with aggregate formation.

B. Effective pair potentials

  • In Figure 6, the authors present the effective potentials obtained from the site-site pair distribution functions.
  • Note that even if in gBB long range correlations due to nanostructure organization are clearly not visible, there are long range repulsions in the BB effective potential, which are reflected in the prepeak in SBB as an indication of IRO.
  • Note that given the large size of the B-sites, the authors have retained the repulsive part of the bare LJ interaction in order to account for the repulsive component of the effective potential.
  • The parameters of the fit are collected in Table III.
  • The range parameter a3 grows considerably with the density, reflecting the increase of intermediate range ordering as the total density is increased.

C. Cluster analysis

  • In order to go beyond the mere qualitative information provided by simulation snapshots and the two-body level information furnished by pair distribution functions or sitesite structure factors, the authors have also performed a geometric cluster analysis on the B sites and the C monomers, using different values for the link distance rcl.
  • The authors will use various values of rcl in the range 10-12 Å for B-sites and C monomers and 6-8 Å for A-sites.
  • Specifically, the authors have calculated the normalized cluster size distribution, N(s), as proposed by Stauffer.45.
  • This is a typical behavior of a non-associating fluid, in which instantaneous clusters are created and destroyed as particles explore their configurational space.

D. The effect of charges

  • The authors previous results have shown that microheterogeneity or stable intermediate range order can be induced by competing short range interactions in a simple mixture model of dimers and monomers.
  • In contrast, Bsites form a percolating bicontinuous structure coexisting with Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions.
  • In spite of the fact that these two effective interactions result from the coarse graining of many body effects, the dominant role of electrostatic interactions already reflected in the partial structure factors leads to surprisingly similar effective potentials when charges are present.
  • (a) Charge dependence of the partial structure factors for A (top), B , and C particles.

IV. CONCLUSIONS

  • This in turn translates into the characteristic presence of a prepeak in the site-site structure factors.
  • These features are found both in simulation and in the integral equation results.
  • The effective site-site potentials extracted from the pair distribution functions by means of an IMC and integral equation approach display the characteristic features of the SALR interactions, with the repulsive long range increasing as the total density (and hence the ordering) increases.
  • It is interesting to note that the appearance of a prepeak in the wide angle scattering experiments and computer simulations of RTILS have been a subject of much investigations47,48 and have been related to the charge ordering and the subsequent appearance of segregated charged and uncharged molecular domains.
  • On the other hand, their simple model when reduced to two dimensions most likely will also give rise to more complex structures, which in three dimensions are hindered by entropic effects.

Did you find this useful? Give us your feedback

Figures (13)

Content maybe subject to copyright    Report

THE JOURNAL OF CHEMICAL PHYSICS 143, 084501 (2015)
Pattern formation in binary fluid mixtures induced by short-range
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, E-28006 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 short-range potentials show
the existence of pattern formation (microheterogeneity), mostly due to depletion forces away from
the demixing region. Eective site-site 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 short-range attractive and a long-range 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 find spherical droplets, sheets, or
tubes. In some cases, the patterns appear as transient states
due to energy or mass fluctuations 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 self-organizing capability of
soft and fluid matter is critical for a wide panoply of appli-
cations of great relevance nowadays. These self-assembly
mechanisms play a crucial role in processes involving protein
solutions in food products,
4,5
therapeutic monoclonal anti-
bodies,
68
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 fluid,
11
a system noto-
rious for its theoretical interest. On the other hand, the addi-
tion of non-adsorbing polymers to the colloidal solution typi-
cally activates an attractive inter-particle interaction, due to the
depletion mechanism. Moreover, changing the concentration
and molecular weight of the polymer, the attraction range
and strength of the colloid-colloid 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 meta-stable
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 colloid-polymer
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 Derjaguin-Landau-Verwey-Overbeek (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
Lennard-Jones (LJ) plus a Yukawa interaction
2,22,23
in order to
model the spontaneous emergence of microstructured patterns
in fluids. 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 air-water interface.
This potential has been studied in depth in model systems, both
in bulk and in confinement,
2529
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
isothermal-isobaric ensembles. We will address the emergence
of intermediate range order (IRO)analyzing the behaviorof the
0021-9606/2015/143(8)/084501/11/$30.00 143, 084501-1 © 2015 AIP Publishing LLC
Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49

084501-2 Bores et al. J. Chem. Phys. 143, 084501 (2015)
partial and concentration-concentration 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 eective interactions
from the pair correlation functions of the simulated mixtures.
For comparison, another set of eective 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 eect leading to
the pattern formation (microheterogeneity, or microstructure
segregation at the nanoscale) translates into the appearance of
eective 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 eective 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 site-site 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 briefly summarize the
methodology. In Section III, we introduce our most significant
results. Conclusions and future prospects are to be found in
Section IV.
II. MODEL AND METHODS
Our model consists in an equimolar fluid mixture of two
dierent species, a two-site dimer AB and a monomer C. The
dimers are represented by a two center LJ site-site 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 site-site 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 coarse-grained 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 nano-structural 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 non-polar tail, by the larger site B. This implies that AA and
CC interactions will be repulsive, BB and AC are attractive,
finally 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. Lennard-Jones 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 fixed 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 non-polar 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 defining 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 eect
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,
3436
in the canonical and isothermal-isobaric ensem-
bles using a Nose-Hoover thermostat and barostat.
37
Our sam-
ples contained 16 384 particles (samples of up to 65 536 parti-
cles were investigated and no significant 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 vapor-liquid 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
isothermal-isobaric simulations and analyzed the volume fluc-
tuation of the samples. In this way, one can avoid those states
that lie inside the liquid-vapor spinodal. Moreover, one can
compute the partial structure factors, defined 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 dierent particles and g
i j
is the atom-atom
Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49

084501-3 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 k-sampling.
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 concentration-concentration structure factor introduced by
Bathia and Thornton,
38
for which we have defined
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 dierent
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 vapor-liquid transition, one can ana-
lyze the corresponding k-dependent linear response suscepti-
bility in density fluctuations, 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 defined 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 vapor-liquid
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
two-phase 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 final 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 site-site eective 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 eective 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
dierent, but in the case of emergence of intermediate range,
order should exhibit some common features.
C. RISM integral equation
The site-site correlations are obtained by solving the usual
set of 2 equations, the site-site Ornstein-Zernike (SSOZ) equa-
tion and the closure equation, which we choose here to be the
site-site hypernetted (SS-HNC) 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 site-site 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 SS-HNC 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/rights-and-permissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49

084501-4 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 defined
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 r-grid 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 eective potentials which
would correspond to the equivalent one-component represen-
tation of the system. This is achieved by imposing the pair
correlation function to be the desired site-site correlation,
namely, g(r) = g
X X
(r), in the set of the two integral equations
for the 1-component system and solve these equations for the
direct correlation function and eective pair interaction. The
direct correlation function can be obtained through the OZ
equation for 1-component 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 eective 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 eective
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 eect of the molecular geom-
etry on the nanostructure formation changing the diameter of
σ
BB
. We have first considered σ
BB
= 8 Å, 9 Å, 10 Å, and 12 Å.
Some snapshots of configurations 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 first peak
of S
BB
(k) exceeds 2.7, which according to the Hansen-Verlet
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 eect, since these are
restricted to occupy the holes between the large B particles.
These eects 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 B-sites 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 first 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 first on the g
AA
pair
distribution function, one first appreciates the large exclusion
hole after the first layer, which is a simple consequence of the
large size of B-sites. Obviously, the exclusion hole grows with
the size of the B-sites, as can be seen when comparing figures
FIG. 1. Snapshots of configurations for total particle density ρ = 0.001 25 Å
3
and temperature T = 226.45 K for two B-site diameters. As the size of
B-sites grows, C monomers cluster in the cavities formed by the B-sites due to excluded volume eects. All other diameters and total density are kept fixed,
σ
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/rights-and-permissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49

084501-5 Bores et al. J. Chem. Phys. 143, 084501 (2015)
FIG. 2. The figures 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 dashed-dotted curves corre-
spond to integral equation calculations.
on the left and right columns. Correlations between A-sites
extend up to five σ
AA
, and the width of the g
CC
correlation
is 2σ
CC
. These features hint at the presence of some degree
of IRO. B-B correlations (graphs in the middle row) behave
like those of a dense fluid, and no apparent sign of clustering
or IRO is evident. In contrast, the wide first peak of g
CC
is
characteristic of clusters of particles confined in cavities, in
this case formed by B-sites. This eect, as mentioned before,
is maximized for the largest σ
BB
. We will see later that these
clusters of partly occluded C-particles 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 B-sites. 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 B-site structure factor shifts to lower k-values
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 concentration-concentration
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
fluctuations inducing spatial inhomogeneities, but no demix-
ing transition. In Figure 5, we have plotted the k-dependent
susceptibility corresponding to density fluctuations. 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 dierent. As density is
decreased, the susceptibility (i.e., the isothermal compress-
ibility) grows, an indication of the vicinity of a vapor-liquid
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 two-phase region. Our
analysis indicates that the thermodynamic conditions we have
simulated can be considered equilibrium states. Moreover, we
have confirmed that the results do not have a significant sample
size dependence, by which metastability can also be ruled out.
The site-site correlation functions and structure factors
obtained from the RISM theory are represented in dashed
lines in Figs. 2-3. 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 significant dierences are seen for the structure factors
in Fig. 3. Integral equations tend to exaggerate concentration
fluctuations 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 fluctuations 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/rights-and-permissions. Downloaded to IP: 161.111.20.160 On: Thu, 31 Mar
2016 12:32:49

Citations
More filters
Journal ArticleDOI
TL;DR: It is shown that the presence, or absence, of a radiation scattering pre-peak is principally related to the symmetry breaking, or not, of the global charge order, induced by the peculiarities of the molecular shapes.
Abstract: The structural properties of ionic liquids and alcohols are viewed under the charge ordering process as a common basis to explain the peculiarity of their radiation scattering properties, namely the presence, or absence, of a scattering pre-peak. Through the analysis of models, it is shown that the presence, or absence, of a radiation scattering pre-peak is principally related to the symmetry breaking, or not, of the global charge order, induced by the peculiarities of the molecular shapes. This symmetry breaking is achieved, in practice, by the emergence of specific types of clusters, which manifests how the global charge order has changed into a local form. Various atom–atom correlations witness the symmetry breaking induced by this re organization, and this is manifested into a pre-peak in the structure factor. This approach explains why associated liquids such as water do not show a scattering pre-peak. It also explains under which conditions core-soft models can mimic associating liquids.

35 citations

Journal ArticleDOI
TL;DR: In this paper, the authors perform an extensive computational study of binary mixtures of water and short-chain alcohols resorting to two-scale potential models to account for the singularities of hydrogen bonded liquids.
Abstract: We perform an extensive computational study of binary mixtures of water and short-chain alcohols resorting to two-scale potential models to account for the singularities of hydrogen bonded liquids. Water molecules are represented by a well studied core softened potential which is known to qualitatively account for a large number of water’s characteristic anomalies. Along the same lines, alcohol molecules are idealized by dimers in which the hydroxyl groups interact with each other and with water with a core softened potential as well. Interactions involving non-polar groups are all deemed purely repulsive. We find that the qualitative behavior of excess properties (excess volume, enthalpy, and constant pressure heat capacity) agrees with that found experimentally for alcohols such as t-butanol in water. Moreover, we observe that our simple solute under certain conditions acts as a “structure-maker,” in the sense that the temperature of maximum density of the bulk water model increases as the solute is add...

20 citations

Journal ArticleDOI
TL;DR: It is proved that in a suitable regime minimizers are periodic stripes, in any space dimension.
Abstract: In this paper we study pattern formation for a physical local/nonlocal interaction functional where the local attractive term is given by the $1$-perimeter and the nonlocal repulsive term is the Yukawa (or screened Coulomb) potential. This model is physically interesting as it is the $\Gamma$-limit of a double Yukawa model used to explain and simulate pattern formation in colloidal systems \cite{BBCH,CCA,IR,GCLW}. Following a strategy introduced in~\cite{DR} we prove that in a suitable regime minimizers are periodic stripes, in any space dimension.

13 citations

Journal ArticleDOI
TL;DR: In this paper, the authors studied the ideality of methanol-ethanol mixtures under ambient conditions of temperature and pressure, and revealed two types of ideality, one related to simple disorder, and another found in complex disorder mixtures of associated liquids.
Abstract: Methanol-ethanol mixtures under ambient conditions of temperature and pressure are studied by computer simulations, with the aim to sort out how the ideality of this type of mixtures differs from that of a textbook example of an ideal mixture. This study reveals two types of ideality, one which is related to simple disorder, such as in benzene-cyclohexane mixtures, and another found in complex disorder mixtures of associated liquids. It underlines the importance of distinguishing between concentration fluctuations, which are shared by both types of systems, and the structural heterogeneity, which characterises the second class of disorder. Methanol-1propanol mixtures are equally studied and show a quasi-ideality with many respect comparable to that of the methanol-ethanol mixtures, hinting at the existence of a super-ideality in neat mono-ol binary mixtures, driven essentially by the strong hydrogen bonding and underlying hydroxyl group clustering.

10 citations

Journal ArticleDOI
TL;DR: Binary mixtures of hard spheres with different diameters and square-well attraction between different particles are studied by theory and Monte Carlo simulations in this article, where the mesoscopic theory, local fl...
Abstract: Binary mixtures of hard spheres with different diameters and square-well attraction between different particles are studied by theory and Monte Carlo simulations. In our mesoscopic theory, local fl...

4 citations

References
More filters
Journal ArticleDOI
TL;DR: In this article, three parallel algorithms for classical molecular dynamics are presented, which can be implemented on any distributed-memory parallel machine which allows for message-passing of data between independently executing processors.

32,670 citations

Journal ArticleDOI
TL;DR: In this paper, the Fourier transform of the pair correlation function is used to calculate the structure factor of a reference system in which the intermolecular forces are entirely repulsive and identical to the repulsive forces in a Lennard-Jones fluid.
Abstract: The different roles the attractive and repulsive forces play in forming the equilibrium structure of a Lennard‐Jones liquid are discussed. It is found that the effects of these forces are most easily separated by considering the structure factor (or equivalently, the Fourier transform of the pair‐correlation function) rather than the pair‐correlation function itself. At intermediate and large wave vectors, the repulsive forces dominate the quantitative behavior of the liquid structure factor. The attractions are manifested primarily in the small wave vector part of the structure factor; but this effect decreases as the density increases and is almost negligible at reduced densities higher than 0.65. These conclusions are established by considering the structure factor of a hypothetical reference system in which the intermolecular forces are entirely repulsive and identical to the repulsive forces in a Lennard‐Jones fluid. This reference system structure factor is calculated with the aid of a simple but accurate approximation described herein. The conclusions lead to a very simple prescription for calculating the radial distribution function of dense liquids which is more accurate than that obtained by any previously reported theory. The thermodynamic ramifications of the conclusions are presented in the form of calculations of the free energy, the internal energy (from the energy equation), and the pressure (from the virial equation). The implications of our conclusions to perturbation theories for liquids and to the interpretation of x‐ray scattering experiments are discussed.

4,462 citations

Journal ArticleDOI
TL;DR: In this article, the scaling theory of phase transition has been used to explain percolation through the cluster properties; it can also be used as an introduction to critical phenomena at other phase transitions for readers not familiar with scaling theory.

1,763 citations

Journal ArticleDOI
TL;DR: As the length of the alkyl chain increases, the nonpolar domains become larger and more connected and cause swelling of the ionic network, in a manner analogous to systems exhibiting microphase separation.
Abstract: Nanometer-scale structuring in room-temperature ionic liquids is observed using molecular simulation. The ionic liquids studied belong to the 1-alkyl-3-methylimidazolium family with hexafluorophosphate or with bis(trifluoromethanesulfonyl)amide as the anions, [Cnmim][PF6] or [Cnmim][(CF3SO2)2N], respectively. They were represented, for the first time in a simulation study focusing on long-range structures, by an all-atom force field of the AMBER/OPLS_AA family containing parameters developed specifically for these compounds. For ionic liquids with alkyl side chains longer than or equal to C4, aggregation of the alkyl chains in nonpolar domains is observed. These domains permeate a tridimensional network of ionic channels formed by anions and by the imidazolium rings of the cations. The nanostructures can be visualized in a conspicuous way simply by color coding the two types of domains (in this work, we chose red = polar and green = nonpolar). As the length of the alkyl chain increases, the nonpolar domai...

1,668 citations

Journal ArticleDOI
TL;DR: In this paper, Monte Carlo computations have been performed in order to determine the phase transitions of a system of particles interacting through a Lennard-Jones potential, and an indirect determination of the phase transition of the hard-sphere gas is made which is essentially in agreement with the results of more direct calculations.
Abstract: Monte Carlo computations have been performed in order to determine the phase transitions of a system of particles interacting through a Lennard-Jones potential. The fluid-solid transition has been investigated using a method recently introduced by Hoover and Ree. For the liquid-gas transition a method has been devised which forces the system to remain always homogeneous. A comparison is made with experiment in the case of argon. An indirect determination of the phase transition of the hard-sphere gas is made which is essentially in agreement with the results of the more direct calculations.

1,101 citations

Frequently Asked Questions (2)
Q1. What are the contributions mentioned in the paper "Pattern formation in binary fluid mixtures induced by short-range competing interactions" ?

Bores et al. this paper showed that a simple mixture of AB dimers and C monomers, with short range attractive and repulsive interactions designed so as to mimic the interactions present in RTILs, can give rise to the presence of nanostructural order in the form of micro-segregation in bicontinuous structures. 

This is certainly a problem relevant to the behavior at interfaces which the authors intend to address in the future.