Liquid-liquid equilibria for soft-repulsive particles: improved equation of state and methodology for representing molecules of different sizes and chemistry in dissipative particle dynamics.

Liquid–liquid equilibria for soft-repulsive particles: Improved equation of state and methodology for representing molecules of different sizes and chemistry in dissipative particle dynamics Thilanga P. Liyana-Arachchi, Sumanth N. Jamadagni, David Eike, Peter H. Koenig, and J. Ilja Siepmann Citation: The Journal of Chemical Physics 142, 044902 (2015); doi: 10.1063/1.4905918 View online: http://dx.doi.org/10.1063/1.4905918 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/142/4?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Mixtures of protic ionic liquids and molecular cosolvents: A molecular dynamics simulation J. Chem. Phys. 140, 214502 (2014); 10.1063/1.4879660 Liquid-vapor phase diagram and surface properties in oppositely charged colloids represented by a mixture of attractive and repulsive Yukawa potentials J. Chem. Phys. 138, 054507 (2013); 10.1063/1.4789915 Phase separation in H 2 O : N 2 mixture: Molecular dynamics simulations using atomistic force fields J. Chem. Phys. 126, 044510 (2007); 10.1063/1.2431171 Isomolar semigrand ensemble molecular dynamics: Development and application to liquid-liquid equilibria J. Chem. Phys. 122, 054504 (2005); 10.1063/1.1839172 Atomistic molecular dynamics simulation of diffusion in binary liquid n-alkane mixtures J. Chem. Phys. 116, 7656 (2002); 10.1063/1.1466472

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 164.15.128.33 On: Tue, 11 Aug 2015 04:07:27

THE JOURNAL OF CHEMICAL PHYSICS 142, 044902 (2015)

Liquid–liquid equilibria for soft-repulsive particles: Improved equation of state and methodology for representing molecules of different sizes and chemistry in dissipative particle dynamics Thilanga P. Liyana-Arachchi,1,a) Sumanth N. Jamadagni,2,a),b) David Eike,2 Peter H. Koenig,2 and J. Ilja Siepmann1,3,b) 1

Department of Chemistry and Chemical Theory Center, University of Minnesota, 207 Pleasant Street SE, Minneapolis, Minnesota 55455, USA 2 Computational Chemistry, Modeling and Simulation, The Procter & Gamble Company, 8611 Beckett Road, West Chester, Ohio 45069, USA 3 Department of Chemical Engineering and Materials Science, University of Minnesota, 421 Washington Avenue SE, Minneapolis, Minnesota 55455, USA

(Received 16 October 2014; accepted 2 January 2015; published online 28 January 2015) Three developments are presented that significantly expand the applicability of dissipative particle dynamics (DPD) simulations for symmetric and non-symmetric mixtures, where the former contain particles with equal repulsive parameter for self-interactions but a different repulsive parameter for cross-interactions, and the latter contain particles with different repulsive parameters also for the self-interactions. Monte Carlo and molecular dynamics simulations for unary phases covering a wide range of repulsive parameters and of densities for single-bead DPD particles point to deficiencies of the Groot and Warren equation of state (GW-EOS) [J. Chem. Phys. 107, 4423 (1997)]. A revised version, called rGW-EOS, is proposed here that is significantly more accurate over a wider range of parameters/densities. The second development is the generalization of the relationship between the Flory-Huggins χ parameter and the repulsive cross-interaction parameter when the two particles involved have different molecular volumes. The third aspect is an investigation of Gibbs ensemble Monte Carlo simulation protocols, which demonstrates the importance of volume fluctuations and excess volumes of mixing even for equimolar symmetric mixtures of DPD particles. As an illustrative example, the novel DPD methodology is applied to the prediction of the liquid–liquid equilibria for acetic anhydride/(n-hexane or n-octane) binary mixtures. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4905918]

I. INTRODUCTION

The dissipative particle dynamics (DPD) methodology was introduced by Hoogerbrugge and Koelman to investigate mesoscale phenomena of complex fluids.1,2 In addition to purely repulsive, conservative forces between DPD particles, such dynamical simulations also involve dissipative and random forces acting on the particles to ensure the correct hydrodynamic behavior.3 Building on this, Groot and Warren (GW) developed a rigorous framework for DPD as a mesoscopic simulation tool for complex chemical systems by establishing an equation of state and a link between DPD crossinteraction parameters and Flory-Huggins χ parameters.4 Enabled by the GW framework, DPD simulations have been utilized to explore homopolymer and blockpolymer melts, surfactants, and a wide range of self-assembly and phase separation applications.5–19,21–28 Building upon the DPD framework, Backer et al.,20 Pivkin and Karniadakis,21 and Spaeth et al.24 have described approaches for systematically changing the coarse-graining degree and for enabling multia)T. P. Liyana-Arachchi and S. N. Jamadagni contributed equally to this

work.

b)Authors to whom correspondence should be addressed. Electronic ad-

dresses: [email protected] and [email protected] 0021-9606/2015/142(4)/044902/13/$30.00

scale simulations, and Pagonabarraga and Frenkel29 and Warren30 introduced many-body extensions to allow for vapor–liquid coexistence. The soft repulsive DPD potential allows for rapid equilibration and exploration of phase space and offers an attractive alternative to traditional allatom force fields with their steeply varying Lennard-Jones or Buckingham repulsive-attractive potentials. The soft-core DPD potentials have also found use in Monte Carlo (MC) simulations.5,6,10,14,24,25 Here, it should be noted that similar gains in efficiency can also be realized by discontinuous molecular dynamics (DMD) simulations where discontinuous potentials allow for the use of collision-based propagation algorithms.31–33 To significantly expand the applicability of the DPD methodology for real-world applications, it is essential to formulate a systematic way to parametrize the soft-core potentials between DPD particles for multicomponent mixtures that accurately reflects both the size and chemical differences between the molecules or molecular fragments that the DPD particles represent. Capturing size differences between particles is important for systems where (entropic) steric effects contribute to aggregate size and shape distributions, e.g., the delicate interplay of surfactant head and tail group sizes for the formation of micelles and bilayers. The ability to parametrize DPD

142, 044902-1

© 2015 AIP Publishing LLC


044902-2

Liyana-Arachchi et al.

J. Chem. Phys. 142, 044902 (2015)

interaction parameters completely in-silico, without relying on experimental data, would be extremely advantageous. For coarse-grained models, there exist no simple mixing rules for determining cross-interaction parameters, such as the widely used Lorentz-Berthelot rules for Lennard-Jones potentials. For a system consisting of M types of DPD beads, O(M 2/2) interactions need to be parametrized. This requires a simplified and fast approach for parameterizing interactions between different particle types. We present such a methodology in this paper.

II. BACKGROUND

In DPD simulations, the particles interact via a softrepulsive, conservative potential ( )  aij r ij 2    1 − if 0 ≤ r ij ≤ r cut , (1) u(r ij) =  2 r cut    0 if r > r ij cut  where r ij is the distance between two particles of types i and j, r cut is the potential truncation distance, and the aij parameter is a positive number (in units of energy) that governs the strength of the repulsive interaction for particles of types i and j. In dimensionless DPD units, r cut is usually set to 1.00 for all particles, but when representing real chemical systems, it is often more convenient to relate r cut to a reference volume (see Sec. IV A). In a molecular dynamics implementation of the DPD methodology, the particles also interact with additional dissipative and random forces to ensure that the temperature is maintained and the fluctuation-dissipation theorem is satisfied.3,4 For unary, single-bead DPD fluids interacting via Eq. (1), Groot and Warren proposed the following equation of state (GW-EOS) by empirically fitting to simulation data:4 PGW(ρ, T) = ρkBT + αaii ρ2,

(2)

where P is the pressure, ρ is the particle number density given as the number of DPD particles in a unit volume (r cut3), T is the temperature, and kB is Boltzmann’s constant. In dimensionless DPD units, kBT is usually set to 1.00 for all systems. Groot and Warren4 obtained α = 0.101 ± 0.001 for the excess pressure arising from the repulsive, pair-wise additive interactions based on three sets of simulations spanning 1.00 ≤ ρ ≤ 8.00 for aii = 15.0, 1.50 ≤ ρ ≤ 8.00 for aii = 25.0, and 2.00 ≤ ρ ≤ 8.00 for aii = 30.0. Using the relationship between P(ρ, T), u(r), and the radial pair distribution function, g(r), for a homogeneous isotropic fluid,34 ( )  ρ2 ∞ ∂u P(ρ, T) = ρkBT − rg(r) 4πr 2dr 6 0 ∂r ( )  rc ( ) 2π r = ρkBT + ρ2 r 3g(r) aii 1 − dr. (3) 3 0 rc  r =1 Maiti and McGrother35 found that α = (2π/3) 0 c r 3g(r)(1 −r/r c) is only a weak function of aii and ρ and is approximately equal to 0.101 over the limited range of simulation conditions assessed. For the case of chemically dissimilar particles of the same molecular volume, Groot and Warren4 showed that a linear relationship holds between aij for the repulsive crossinteraction and the Flory-Huggins χ parameter. Building on

this, Maiti and McGrother35 showed that the DPD interaction parameters could be related to the Hildebrand solubility parameters and that properly parametrized DPD simulations yield the correct interfacial tension for different levels of coarse graining. Along similar lines, Travis and coworkers36 used the cohesive energy density to determine aii and the correspondence between DPD and regular solution theory instead of the χ parameter to determine aij. However, in all these applications of the DPD methodology to (usually binary) mixtures, the molecular volume of different particle types was not considered a target property for parameterization. We demonstrate how to overcome the limitation to equal-sized particles in Secs. IV and V. Specifically, we derive and validate a revised GW-EOS, called the rGW-EOS, with improved accuracy for a larger range of densities and aii parameters. We also show how to parametrize the cross-interaction term for athermal mixtures of DPD particles of different molecular volumes. We then extend the linear relationship between the Flory-Huggins χ parameter and the aij cross-interaction parameter to the general case when the interacting particles have different molecular volumes.

III. SIMULATION METHODS

To generate the large sets of simulation data for single-bead DPD particles presented in this manuscript, we relied predominantly on MC simulations in various ensembles, but molecular dynamics simulations were also utilized for a more limited set of conditions (see supplementary material54). The in-house MCCCS-MN (Monte Carlo for Complex Chemical SystemsMinnesota) software was used for the DPD simulations.37 To investigate the limitations of the GW-EOS and aid in the development of the rGW-EOS, Metropolis MC simulations in the canonical ensemble38,39 were used, and single-phase MC simulations in the canonical and isobaric-isothermal ensembles40 were used to probe athermal mixtures. For constantNVT simulations, only center-of-mass translational moves are needed, whereas both center-of-mass translations and volume displacements are employed in constant-N PT simulations. For all simulations, the maximum displacements for translational and volume moves were adjusted during the equilibration period to yield acceptance rates close to 50% and 40%, respectively. The total number of DPD particles was set to Ntotal = 500 and 1000 for unary systems and binary mixtures, respectively. The fraction of volume moves was adjusted to yield one accepted move per MC cycle, where a cycle consists of Ntotal moves. Four independent simulations were carried out at each state point with equilibration and production periods consisting of at least 105 and 2 × 105 MC cycles, respectively. To compute the infinite-dilution excess chemical potential or the Gibbs free energy of transfer for a particle of type i from a reference vapor phase to a liquid phase consisting of particles j, we performed MC simulations in a pseudo-osmotic variant of the Gibbs ensemble,41 where Ni = 1 and N j = 500, volume moves were applied only to the liquid phase and swap moves were only performed for the i-type particle. For the latter move type, the configurational-bias MC algorithm was used.42,43 Again, the fraction of moves was adjusted to yield


044902-3


J. Chem. Phys. 142, 044902 (2015)

approximately one accepted volume and one accepted swap move per MC cycle. Four independent simulations were conducted. The Gibbs free energy of transfer was calculated according to44,45 ρi ex * liq, j + , ∆Gtrans = µ = −RT ln i, j i, j i , ρvap -

(4)

where ρiliq, j and ρivap are the average number densities of the i-type particle in the liquid and vapor phases, respectively, and R is the molar gas constant. Liquid–liquid coexistence data for binary mixtures were determined using both constant-volume (NVT) and constantpressure (N PT) Gibbs ensemble Monte Carlo (GEMC) simulations.41,46 Two simulation boxes in thermodynamic contact but not sharing explicit interfaces were used for the GEMC simulations with up to three types of MC moves being employed: center-of-mass translations, volume exchanges (between the two boxes for constant-volume GEMC and with a pressure reservoir for constant-pressure GEMC), and configurational-bias particle swaps between the two simulation boxes. All binary systems contained Ntotal = 1000 DPD particles. At least four independent simulations were carried out at each state point with equilibration and production periods consisting of at least 2 × 105 and 5 × 105 MC cycles, respectively. For the calculation of phase equilibria for DPD systems near the critical region, GEMC simulations are better suited than molecular dynamics simulations with a single elongated simulation box containing two phases (and two interfaces).5

IV. RESULTS AND DISCUSSION A. Representing molecules or molecular fragments of different sizes in multicomponent DPD simulations

By size, we refer to the molecular volume of a singlebead DPD particle (i.e., volume in units of Å3 occupied on average by a single molecule or a group of specific molecules) or the corresponding fragment volume of a specific fragment belonging to a molecule that is represented by multiple DPD beads (e.g., a three or four carbon segment in a long chain alkane or a polydimethylsiloxane monomer in a silicone polymer may constitute a fragment). This manuscript is focused on the former case, i.e., single-bead DPD particles. In contrast to a Lennard-Jones particle characterized by a size parameter, σ, and a well depth, ϵ, there is no explicit “size” of a particle described by the one-parameter DPD potential. However, we describe below how the pure-phase particle density, ρ0i , and the DPD repulsive parameter, aii, can be suitably chosen to represent molecules (or groups of molecules or fragments of molecules) with different sizes. Here, we should note that the issue of size has previously been addressed in the context of the degree of coarse graining, and methodologies have been developed that involve both changes in the cutoff radius and repulsive parameter to change the length scale of the system.20,21,24 As is commonly done in DPD simulations, we often express physical quantities such as density, temperature, pressure, and energy in dimensionless units. To map DPD particles

onto actual molecules, we choose a reference particle with a known volume. Water, being the most common solvent, is a good choice but, if one is not interested in aqueous systems, then another molecule (or the geometric mean for the molecules in the mixture) would be a suitable alternative. To increase the level of coarse graining, it is common to group 3-4 water molecules into a single water-like DPD particle.9,23,26,27 This makes the molecular volume of the water-like particle comparable to that of other solvent molecules (e.g., ethanol or acetonitrile) or of common molecular fragments. The molecular volume of a single water molecule at ambient conditions is ≈30.0 Å3, giving vref = 3×30.0 = 90.0 Å3 (choosing the volume of such a water-like particle as the reference volume). We also have to choose a reference number density, ρref , that determines r cut in physical units and allows one to set the coarse-graining length scale for the DPD simulations. As will be seen later, ρref needs to be selected carefully so that a wide range of different molecular or fragment volumes can be considered. In this work, we chose ρref = 5.00 for most applications which falls within the density range of 3.00 to 6.00 recommended by Groot and Warren.4 For our choice of vref = 90.0 Å3 and ρref = 5.00, one finds r cut = (vref ρref )1/3 = 7.66 Å. Groot and Warren4 also suggested the following relationship between the repulsive aii parameter and the DPD particle number density by matching the isothermal compressibility as calculated from the GW-EOS to that of the water-like particle at ambient conditions aii ρ κ −1 − 1 = ≈ 75. (5) kBT 2α Thus, for the reference system (with kBT set to 1), we obtain aref = 75/ρref = 15.0 and inserting these values into the GWEOS (Eq. (2)) yields a predicted pressure of PGW = 42.9. This equation sets the degree of coarse-graining to unity21 which is appropriate for the present work because the focus is on representing each molecule by a single-bead DPD particle (irrespective of how many interaction sites such a molecule would possess in a molecular mechanics force field). Following Kacar et al.,25 we set the pure-phase number density of each bead type in a multicomponent mixture as follows: vref ρref ρ0i = . (6) vi Accounting for these size differences is important to capture packing or steric effects in simulations (for example, straight versus branched chain surfactants in a bilayer). With this approach, particles that represent molecules or molecular fragments larger than the reference particle will have a lower pure-phase number density and vice versa. For example, acetic anhydride has a molecular volume of 157 Å3. If it is represented using a single particle, with the water-like particle as the reference particle, the acetic anhydride particle density would be 2.87 (=90.0×5.00/157). Molecular volumes can be accurately calculated from macroscopic liquid densities, which can be predicted from a number of computational tools, and Eq. (6) provides a facile route to incorporate that information into a DPD simulation. With ρref chosen, we set aref using Eq. (5) and the reference bead’s compressibility, and use the GW-EOS (or the


044902-4


J. Chem. Phys. 142, 044902 (2015)

rGW-EOS described later) to calculate the reference pressure, PEOS. If different DPD particles occurring in a mixture have different pure-phase number densities, then their repulsive selfinteraction parameters (aii) must now be suitably adjusted to yield the same pressure predicted by the EOS for the reference system. For the GW-EOS, this common-pressure constraint yields aii =

PGW − ρ0i kBT α(ρ0i )2

.

(7)

Equations (6) and (7) together provide a self-consistent method to parametrize different molecules and molecular fragments as particles in DPD simulations of unary systems and multi-component mixtures, but it inherently assumes that the compressibility of the reference system also applies to all other systems. B. Limitations of the GW-EOS and development of the rGW-EOS

To check the accuracy of the GW-EOS over a wide range of conditions, we performed simulations of single-component, single-bead DPD fluids covering a range of molecular volumes (i.e., different particle number densities) with aii computed from Eq. (7). To obtain PGW, reference densities of 3.00 and 5.00 were chosen for this part; Eqs. (5) and (2) yield aref = 25.0 and 15.0 and PGW = 25.7 and 42.9 for these systems, respectively. Fig. 1 shows a comparison of the pressures computed from simulations and the target pressure predicted by the GWEOS. While the GW-EOS appears reasonably accurate at high densities, there are rather significant deviations at lower densities as already noticed by Groot and Warren.4 The relative deviations are 3% and 8% for the systems at their ρref of 5.00 and 3.00. Preliminary simulations for liquid–liquid equilibria (LLE) of the soft-core DPD particles that are much more compressible than systems described by Lennard-Jones potentials indicate that such deviations can have a rather significant impact on the phase envelope and, hence, hamper the predictive ability of DPD simulations for particles with different molecular volumes. This observation motivated us to re-visit the

FIG. 1. Comparison of pressures computed from simulations (symbols) to those predicted from the GW-EOS (lines). Data for ρ ref = 5.00 are shown as circles and dashed line and those for ρ ref = 3.00 as squares and dotted line.

development of an EOS for DPD particles and to propose a revised GW-EOS, called the rGW-EOS. As also alluded to by Groot and Warren,4 the ratio of the excess pressure to the squared density should approach the second virial coefficient, B2(aii), in the low-density limit and not the prefactor of αaii found in the GW-EOS. This leads us to consider a more complex EOS that employs two densitydependent switching functions as follows: PrGW(ρ, T) = ρkBT + [ f B2(ρ)B2(aii) + f a (ρ)aii] ρ2,

(8)

where α from the GW-EOS is subsumed into f a (ρ). B2(a) can be computed by numerical integration of the Mayer function34 for the (conservative part of the) DPD potential (Eq. (1))  rc (exp[−u(aii, r)/k BT] − 1)r 2dr. (9) B2(aii) = −2π 0

Here, it was found that the following polynomial fit closely reproduces B2(aii): B2(aii) = 1.705 × 10−8aii5 − 2.585 × 10−6aii4 + 1.556 × 10−4aii3 − 4.912 × 10−3aii2 + 9.755 × 10−2aii.

(10)

To determine suitable functional forms and numerical parameters for f B2(ρ) and f a (ρ), we first performed NVT simulations of DPD fluids (with 500 particles in the simulation box) over a wide range of number densities (from 0.125 to 8.0) and self-interaction parameters (from 10.0 to 40.0). The simulation results are shown in Fig. 2; as suggested by Groot and Warren, 4 the data are presented not as the raw pressure but in the form of the normalized excess pressure, (P− ρk BT)/a ρ2, which accentuates the difference between the different curves. The main part of Fig. 2 illustrates that data for different aii parameters collapse onto each other for ρ ≥ 3.00, but there is a considerable spread at lower densities with larger aii parameters yielding smaller B2 values. Even for ρ ≥ 3.00, the excess pressure continues to increase, albeit very gradually with increasing density (see inset of Fig. 2). The GW-EOS, however, approximates this region as yielding a constant value that is described by α = 0.101 ± 0.001.

FIG. 2. Calculated values of the normalized excess pressure as a function of density for all a ii parameters studied. The inset shows a close-up for the density range commonly considered for DPD simulations with Groot and Warren’s α = 0.101 indicated by the dashed line. The legend provides the a ii values for the different symbols.


044902-5


J. Chem. Phys. 142, 044902 (2015)

To improve the GW-EOS, we have to choose the appropriate functional forms for f B2(ρ) and f a (ρ) in Eq. (8). However, simultaneous optimization of the two switching functions is difficult due to their unknown functional forms and the multiple parameters involved. Instead, we rearrange Eq. (8) to solve for f a (ρ) f a (ρ) =

PrGW − ρk BT − f B2(ρ)B2(aii)ρ2 . aii ρ2

(11)

This provides a route to find an appropriate functional form for f B2(ρ) by choosing f B2(ρ) such that the simulation data for different aii values collapse onto one curve when the right hand side (RHS) of Eq. (11) is plotted as function of ρ (see Fig. 3). Subsequently, one can fit a functional form to the collapsed data to obtain f a (ρ). Qualitatively, f B2(ρ) should approach unity as ρ approaches zero and should decay to zero for moderate to high densities. Empirically, we find the following functional form to be satisfactory for f B2(ρ):

forms, f B2(ρ) decays more rapidly toward zero with increasing pressure than f a (ρ) increases toward its limiting value of c1/c2 = 0.103. Albeit this limiting value is somewhat larger than α in the GW-EOS, it is now closer to the mean-field value of π/30 = 0.1047.20 The actual values of f a at ρ = 5.00 and 7.00 are only 0.0980 and 0.1004, respectively, i.e., they fall below the GW value over the entire range of densities commonly used in DPD simulations. The resulting rGW-EOS for the single-bead DPD fluid is thus given by   B2(aii) aiic1 ρ2 2 + ρ. (14) PrGW(ρ, T) = ρk BT + 1 + ρ3 1 + c2 ρ2

It should be noted that plotting data for the GW-EOS in a similar manner, i.e., using (Psim − ρkBT)/aii ρ2 ≈ 0.101, indicates significant deviations for ρ < 5. The two switching functions are compared in Fig. 3. Due to their different functional

As illustrated in Fig. 4, the rGW-EOS holds very well over most of the aii and ρ space used for the parameterization. At low densities (and correspondingly low pressures), the data for a given ρ and different aii line up as a group and those for ρ = 1.00 show the largest errors (>4% for aii ≥ 35). When grouping data for a given aii, then a characteristic shape emerges for ρ ≤ 0.25, the rGW-EOS very slightly underpredicts Psim; as the density increases to intermediate values, the rGW-EOS overpredicts (with a maximum deviation found at ρ = 1.00 for aii ≥ 20); this is followed by another region where the rGW-EOS slightly underpredicts; finally, the rGWEOS approaches Psim for the higher densities. This rather non-monotonic behavior of the residual illustrates that further gains in accuracy would require a significantly more complex functional form. In contrast, the GW-EOS yields significant deviations over a large fraction of the (aii,ρ) space. The mean unsigned percentage errors (MUPEs) are 11.5% and 0.7% for the GW-EOS and rGW-EOS, respectively, when all data points are considered. Even when only considering the range of densities most commonly used in DPD simulations (3 ≤ ρ ≤ 7),

FIG. 3. The two switching functions of the rGW-EOS as function of density: f B2 (top) and f a (bottom). The symbols in the bottom panel show the simulation data in form of the RHS of Eq. (11). The dashed lines show the corresponding data for the GW-EOS, i.e., 0.101 − f B2 B2(a ii)/a ii.

FIG. 4. Performance of the rGW-EOS and GW-EOS: percentage errors for the rGW-EOS (top) and logarithm of the unsigned percentage error as function of pressure. The different symbols illustrate either the a ii value (top) or the EOS (bottom).

f B2(ρ) =

1 . 1 + ρ3

(12)

As illustrated in Fig. 3, this switching function indeed yields an almost complete collapse of the data for the RHS of Eq. (11) over the full density range. The following expression provides a good fit of these collapsed data: f a (ρ) =

c1 ρ2 1 + c2 ρ2

with

c1 = 0.0802

and

c2 = 0.7787. (13)


044902-6


J. Chem. Phys. 142, 044902 (2015)

then the MUPE for the rGW-EOS remains significantly lower than that for the GW-EOS (0.5% versus 3.7%). C. Limitations of the rGW-EOS

To examine the accuracy of the rGW-EOS, we performed additional MC simulations covering a range of densities for a single-component, single-bead DPD fluid. For this test, we chose a reference system with ρref = 5.00 and aref = 15.0 for which the GW-EOS and rGW-EOS predict pressures, PEOS, of 42.9 and 41.9, respectively. For a desired density, Eq. (7) yields the aii parameter that should yield PEOS for the GWEOS, whereas Eq. (14) needs to be iteratively solved to yield the aii parameter for the rGW-EOS because the second virial coefficient also depends on this parameter. Fig. 5 shows how the actual simulation pressures calculated for the deduced aii parameter deviate from PEOS. As discussed earlier, using the GW-EOS to deduce the self-repulsion parameter results in simulation pressures lower than PEOS, especially at lower densities (for ρ ≤ 3.00, the relative errors are larger than that of 7%). However, since the sign of these deviations for the GW-EOS does not change, simulation of mixtures consisting of particles with different molecular volumes benefit from error compensation. For example, when two immiscible fluids are brought into mechanical contact, the deviations in their densities are smaller than when the corresponding fluids would be probed separately in contact with a pressure bath. Overall, the rGW-EOS performs much better and yields relative errors less than 1% for ρ ≥ 3.3 and a 2% error for ρ = 3.00. The rGW-EOS is least accurate for 2.00 ≤ ρ ≤ 2.75, where the simulation pressures are at least 3% larger than PEOS. At even lower densities, the B2 term in Eq. (14) starts to dominate and the errors decrease again. A possible reason for the reduced accuracy of the rGW-EOS at low densities (ρ ≤ 2.75) is that very large aii parameters are needed to achieve PEOS = 41.9 (see inset of Fig. 5). For example, at ρ = 3.00, the aii parameter is 47.5 for the rGW-EOS (and 43.9 for the GWEOS). This value falls slightly above the maximum value of 40 considered for the development of the rGW-EOS. When

FIG. 5. Percentage errors obtained when using the GW-EOS and rGW-EOS to deduce the self-repulsion parameter needed to match the target pressures of PEOS = 42.9 and 41.9, respectively, calculated with these EOS for a reference system (ρ ref = 5.00 and a ref = 15.0). The inset shows the corresponding a ii parameters.

FIG. 6. Radial distribution functions calculated for DPD fluids at different densities (see legend for line colors) with a ii chosen to yield PEOS = 41.9 obtained from the rGW-EOS for a reference system.

an attempt was made to include data from simulations with aii ≥ 45 in the development of the rGW-EOS, we noticed that the RHS of Eq. (11) did not collapse as well onto one curve (like observed in Fig. 3) for any choice of f B2(ρ) tested. Instead of trying much more complex switching functions or using switching functions that are functions of ρ and a, we studied how the structure of the DPD fluid changes for very large aii values to better understand the physical origin of the limitations of the DPD protocols. Fig. 6 shows the radial distribution functions (RDFs) for DPD fluids at different densities with the aii parameter deduced from the rGW-EOS to match the pressure for the reference fluid (ρref = 5.00 and aii = 15.0). For the large aii parameters needed at low densities (ρ ≤ 3.00), the RDF becomes very structured with three clearly discernible peaks and a large excluded region for r < 0.5. These features indicate that such DPD systems can no longer be characterized as “soft” fluids with significant overlap between particles, but start to resemble more the RDFs common for Lennard-Jones particles. In contrast, the RDFs for the DPD fluids at ρ ≥ 3.50 are much less structured. For example, the heights of the first and second peaks and the depth of the first minimum are only 1.17, 1.03, and 0.93, respectively, for the DPD fluid at ρ = 3.50. In addition, the excluded regions are much smaller and for ρ ≥ 5.00, the values of g(r = 0) are non-zero. Similar features can be observed for the potentials of mean force (PMFs) between particles, w(r) = −k BT ln[g(r)] (see Fig. 7). For ρ ≤ 3.00, the PMFs are non-zero and oscillatory even at separations two times larger than the cutoff (r cut = 1.00) and w(r) > kBT at short separations (r < 0.4). Consequently, the advantages of the soft DPD potential are somewhat diminished when large aii values must be used at low densities to achieve PEOS from a well-behaved reference fluid at ρ = 5.00. When using “larger” particles with lower number densities to represent larger molecules/molecular fragments (i.e., fewer number of interaction sites and greater level of coarse-graining),20,21,24 then smaller time steps have to be used to avoid MD integrator errors due to the stronger forces. The larger PMFs observed for ρ ≤ 3.00 are the cause for the larger deviations in the rGW-EOS (and also in the GWEOS) because of a breakdown of the mean-field approximation


044902-7


J. Chem. Phys. 142, 044902 (2015)

FIG. 7. Potentials of mean force calculated from the RDFs shown in Fig. 6.

underlying the third term of Eq. (14) (and the second term of Eq. (2)). In light of these observations, we recommend that caution be exercised when pushing the coarse-graining level in DPD simulations. For the specific case of a water-like particle at ambient conditions being the DPD reference system with ρref = 5.00 and aref = 15.0, we recommend that molecular volumes of other solvent molecules or molecular fragments that are represented as single DPD bead do not exceed 150 Å3, so that ρ0i > 3.00 (using Eq. (6)). For this density range, the rGWEOS yields relative deviations of less than 2% (see Fig. 5). D. Athermal mixtures

Before investigating mixtures of particles dissimilar in both size and chemistry, we first examined “athermal” mixtures of particles with different molecular volumes where the cross-interaction parameter is suitably chosen.35,36 To reiterate, different “sizes” imply that the beads represent molecules or molecular fragments with different molecular volumes (e.g., ethanol and isopropanol), and consequently, have different pure-phase number densities at the same PEOS. An ideal mixture is characterized by (i) no volume change upon mixing, (ii) no enthalpy change upon mixing, and (iii) an entropy of mixing given only by the combinatorial pre-factor. In contrast, only condition (ii) needs to be satisfied for an athermal mixture. For a binary mixture, the two self-interaction parameters, aii and a j j can be determined from the rGW-EOS with their ρ0i and ρ0j as input that are determined from the molecular volumes. From consideration of the RDFs for binary mixtures (see Eq. (3)), Travis et al.36 deduced an equation for the cross-interaction parameter aij =

aii (ρ0i )2 + a j j (ρ0j )2 2ρ0i ρ0j

.

(15)

By analogy, we use the following ansatz for the rGW-EOS: aij =

h(aii, ρ0i ) aii (ρ0i )2 + h(a j j , ρ0j ) a j j (ρ0j )2 ,  2ρ0i ρ0j h(aii, ρ0i ) h(a j j , ρ0j ) B2(a) f B2(ρ0) h(a, ρ0) = f a (ρ0) + , a

(16)

FIG. 8. Thermodynamic properties for binary mixtures as function of mole 0 > ρ 0 ): fraction of the A-type particle with the smaller molecular volume (ρ A B relative deviation of the calculated pressure from PEOS (top), relative excess internal energy (top middle), volume (bottom middle), and enthalpy (bottom) of mixing. The dashed lines and symbols denote results calculated using Eqs. (15) and (16), respectively. The colors indicate the system specification: 0 = 5.00 and ρ 0 = 3.50, (red circles) ρ 0 (blue squares) ρ ref = ρ A ref = ρ A B 0 0 = 5.00 and ρ B = 3.00, (green up triangles) ρ ref = 5.00, ρ A = 7.07, and 0 = 7.07 and ρ 0 = 3.54. ρ B0 = 3.54, (cyan down triangles) ρ ref = ρ A B

where h(a, ρ0) is the prefactor to the a ρ2 term in the rGW-EOS (Eq. (14)) written in the form P = ρkBT + h(a, ρ0)a ρ2 (with a being either aii or a j j ). To check the validity of Eqs. (15) and (16), we performed NVT and N PT MC simulations for binary mixtures consisting of particles with different molecular volumes (and correspondingly different ρ0i and ρ0j ). Fig. 8 shows the behavior for different ratios of ρ0A/ρ0B and different reference densities. Keeping ρref = 5.00, we explored size ratios of 1.429, 1.6, and 2 and set ρ0A = ρref for the two smaller ratios, but used ρref as the geometric mean of ρ0A and ρ0B for the largest size difference to remain within the recommended density range for DPD simulations (i.e., ρ0B > 3.00). The fourth test case also uses a size ratio of 2, but ρref is now set to ρ0A, the pure-phase density of the smaller particle. For ρref = 5.00 and 7.07 and the usual reference compressibility, the rGW-EOS yields PEOS = 41.9 and 60.4, respectively. For NVT simulations, the relative deviations in the pressure are governed by the deviations of the rGW-EOS for the two pure phases and are nearly linearly dependent on the mixture composition (see top part of Fig. 8). For the mixture with the largest size difference, the relative error in the pressure is smaller for ρref = 5.00, whereas the opposite is true for the relative error in the enthalpy of mixing. The ratio of the average deviation for these two ρref s is somewhat larger for the pressure than for the internal energy and, hence, it appears


044902-8


J. Chem. Phys. 142, 044902 (2015)

that ρref = 5.00 is the slightly better choice. When the data for the excess properties of mixing are compared for different size ratios, then it is clear that larger size ratios lead to larger deviations from ideal mixing. The same observations hold true for N PT simulations. Comparing the data obtained for Eqs. (15) and (16), it appears that using Eq. (15) leads to slightly smaller deviations from ideality for the volume of mixing than using Eq. (16). However, the reverse holds true and to a larger extent for the internal energy and the enthalpy of mixing. Overall, use of Eq. (16) leads to smaller deviations and also holds the advantage of being more consistent with the form of the rGWEOS. Thus, we recommend using Eq. (16) as the starting point to obtain aij for mixtures. However, it should be noted that neither equation truly results in an athermal mixture because they yield a non-zero enthalpy (and internal energy) of mixing. Comparing the relative errors for NVT and N PT simulations, we note that the average deviations are significantly smaller for N PT simulations (e.g., by a factor of about 2 for the enthalpy versus internal energy of mixing). For N PT simulations, a size ratio of 0.6 yields excess properties of mixing with average deviations smaller than 0.2% which is certainly within an acceptable range for DPD simulations. This small error could likely be even further reduced when the densities are selected in a manner that ρref is their geometric mean (i.e., ρ0A = 6.45 and ρ0B = 3.87).

×

E. Cross interactions for thermal mixtures (χ > 0)

In their pioneering work, Groot and Warren also showed that fluid mixtures of DPD particles (interacting via the oneparameter DPD potential and described by the GW-EOS) can be mapped onto Flory-Huggins theory.4 For binary mixtures of DPD particles of the same pure-phase number density, they observed a linear relationship between the Flory-Huggins χ parameter and the repulsive cross-interaction parameter, aij, as follows: χ = cρ [aij − aii],

0 where µex i, j (aij,ρ j ) is the infinite dilution excess chemical potential of particle i in fluid j at density ρ0j with aij describing the cross-interaction (where the self-interaction parameters aii and a j j are determined from the rGW-EOS at a specific PEOS). The task is then to choose aij for a given set of {χ,ρ0i ,ρ0j } such that Eq. (19) is satisfied. 0 To this end, we first calculated µex i, j (aij,ρ j ) over a two0 dimensional grid in terms of ρ j and aij. The pure-phase number density of particle type j was varied from 3.00 to 6.00 in increments of 0.50 where the self-interaction parameter a j j was calculated using the rGW-EOS with ρref = 5.00 and aref = 15.0 resulting in PEOS = 41.9. The cross interaction parameter aij was varied from 10.0 to 90.0 in increments of 5.0. The simulations data are shown in Fig. 9. The excess chemical potential increases monotonically with increasing ai, j and approaches the finite limiting value for transfer of a hard-sphere particle. Increasing ρ0j also leads to an increase in µex i, j because of the increase in the entropic cost for cavity formation. Since using a two-dimensional interpolation for the µex i, j (aij,ρ0j ) data while solving Eq. (19) iteratively for aij can become cumbersome, we searched for an analytical function as replacement. As illustrated in Fig. 9, the following equation yields an accurate fit of the data with a MUPE of 0.3%:   ( ) 0 0 µex i, j (aij,ρ j ) = 1 − exp −aij 0.00323ρ j + 0.00439

(17)

50.9ρ0j 1 + 0.790ρ0j

.

(20)

With this equation, aij can be calculated for any {χ,ρ0i ,ρ0j } by iteratively solving Eq. (19) (we used a simple NewtonRaphson solver). However, within the DPD framework, it may be more advantageous to set aij to a specific value and to use Eq. (20) to calculate the excess chemical potentials that in turn can be used as input for Eq. (19) to calculate χ and then the corresponding temperature from Eq. (18). For ease of use, we also explored whether aij(χ,ρ0i ,ρ0j ) resulting from this numerical procedure can be described by

where the fitting coefficient cρ was found to be 0.286 and 0.689 for ρref = 3.00 and 5.00, respectively.4 Here, we develop a set of equations for the generalized case of mixtures consisting of DPD particles of different sizes (i.e., ρ0i , ρ0j ). We begin by expressing the χ parameter in terms of transfer free energies10 ( ) vref ∆G i, j (T) ∆G j,i (T) + (18) χ≡ 2RT vi vj where ∆G i, j (T) is the (infinite dilution) free energy change for transferring a particle of type i from its pure phase into a pure phase of type j at absolute temperature T, and vi is the molecular volume of the particle of type i. An alternative form of Eq. (18) replaces the transfer free energies with temperaturedependent Hildebrand solubility parameters.28 Equation (18) can be expanded in form of the excess chemical potentials as    1 0 ex 0 ρ0i µex χ= i, j (aij,ρ j ) − µ i,i (aii,ρ i ) 2RT ρref   0 ex 0 + ρ0j µex (19) j,i (aij,ρ i ) − µ j, j (a j j ,ρ j ) ,

FIG. 9. Infinite dilution excess chemical potentials for particles with their cross interaction parameter set to a ij in fluids with different pure-phase densities, ρ 0j at PEOS = 41.9 (and with the corresponding a j j ). The symbols depict the simulation data and the dashed lines are fits according to Eq. (20). See legend for symbol style and color coding.


044902-9


J. Chem. Phys. 142, 044902 (2015)

an analytical function similar to Eq. (17) for the simpler case of equal-sized particles. Replacing the coefficient cρ with a to be determined function, f (χ,ρ0i ,ρ0j ), and replacing interaction parameter aii with aij(χ = 0) from Eq. (16) for the athermal mixture leads to χ = f ( χ,ρ0i ,ρ0j )   h(aii,ρ0i ) aii(ρ0i )2 + h(a j j ,ρ0j ) a j j (ρ0j )2    , × aij −    2ρ0i ρ0j h(aii,ρ0i )h(a j j ,ρ0j )   (21) constraint f (χ,ρ0i ,ρ0j ) = f (χ,ρ0j ,ρ0i ) must specific case of PEOS = 41.9 (from ρref

where the symmetry be satisfied. For the = 5.00 and aref = 15.0), the following equation fits the data obtained from the Newton-Raphson solver very well (see Fig. 10): aij =

h(aii,ρ0i )aii(ρ0i )2 + h(a j j ,ρ0j )a j j (ρ0j )2  2ρ0i ρ0j h(aii,ρ0i )h(a j j ,ρ0j ) +

30.17χ 561.6χ1.667 + 0 0 2.864 . ρ0i ρ0j (ρi ρ j )

(22)

For comparison, if ρ0i = ρ0j = 5.00 and aii = a j j = 15.0, Eq. (22) simplifies to ∆aij = aij − aij(χ = 0) = 1.207χ+0.05568χ1.667 and yields values of 6.849 and 14.65 for χ = 5 and 10, respectively, that are similar to the values of 7.257 and 14.51, respectively, obtained from Eq. (17). For ρ0i = ρ0j = 3.00, the values determined from Eqs. (22) and (17) cannot be compared directly because they correspond to different PEOS states. Liquid–liquid equilibria are known to be quite sensitive to minor variations in the interaction parameters. To probe the effect of using the rGW-EOS with aij parameters determined from Eqs. (19) and (20) versus using the GW-EOS and Eq. (17), we carried out simulations for a system with ρ0i = ρ0j = 5.00. As can be seen from the data presented in Fig. 11, the rGW-EOS yields an upper critical solution temperature (inverse χ parameter) that is about 20% lower than that for the GW-EOS approach because

FIG. 10. Scatter plot showing the quality of the empirical expression (Eq. (22)) to fit ∆a ij(χ, ρ 0i , ρ 0j ).

FIG. 11. Liquid–liquid equilibria for particles with equal molecular volumes (ρ 0i = ρ 0j = 5.00) and like and unlike repulsive interaction parameters determined from the GW-EOS (circles) and rGW-EOS (triangles) approaches. The dashed and dotted lines and filled symbols indicate fits to the scaling law and the corresponding upper critical solution points.

a given χ value corresponds to a smaller ∆aij value for the rGWEOS approach. Correspondingly, the mutual solubilities are found to be larger by 28% and 56% for the rGW-EOS at χ = 5 and 2.75, respectively, than those obtained with the GW-EOS approach. If one routinely wants to carry out DPD simulations for PEOS , 41.9 (i.e., ρref , 5.00 and aref , 15.00), then a similar two-dimensional grid of µex i, j data as shown in Fig. 9 would have to be generated at that specific PEOS and an expression similar to Eq. (22) would need to be parameterized. F. GEMC simulation protocols for liquid–liquid equilibria

Before applying the rGW-EOS and the relationship between the Flory-Huggins χ parameter and the cross-interaction parameter aij to the prediction of the LLE for a specific chemical mixture, we turned our attention to exploring various protocols for GEMC simulations. To this extent, we computed the LLE for a mixture of single-bead, equal-sized DPD particles that was previously investigated by Wijmans et al.10 For this mixture, ρref = ρ0A = ρ0B = 3.00, aAA = aBB = 25.0, and the repulsive cross-interaction parameter, aAB = aAA + ∆a are varied (equivalent to varying χ or the temperature). Wijmans et al.10 suggested that computationally expensive volume moves would not need to be performed for such a symmetric mixture of equal-sized particles when an equimolar system (NA = NB) is used for simulations in the canonical version of the Gibbs ensemble.46 In order to explore this matter and to establish correct GEMC simulation protocols for DPD particles, we performed five sets of GEMC simulations: NVTGEMC simulations without and with volume moves for both equimolar (NA = NB = 500) and non-equimolar (NA = 600 and NB = 400) systems and also N PT-GEMC for an equimolar system. We calculated both the average mole fractions and the mole fraction distributions in the two phases to validate the correct approach for GEMC simulations of DPD particles. For


044902-10


FIG. 12. Liquid–liquid coexistence curves computed using different GEMC protocols for a mixture consisting of two different DPD particles with equal size. The black circles denote the data reported by Wijmans et al.10 The colored symbols and line styles denote simulation data for 1000-particle systems using the following GEMC protocols: N V T ensemble, NA = NB = 500 without volume moves (red squares) and with volume moves (green up triangles); N V T ensemble, NA , NB without volume moves (blue diamonds) and with volume moves (purple down triangles); and N PT ensemble with NA = NB = 500 (magenta stars). The magenta and cyan symbols show data for 4000- and 8000-particle systems, respectively, without volume moves (crosses) and with volume moves (pluses). Uncertainties in the coexistence compositions have been omitted for clarity, but are provided in Table S6.

equimolar systems, Fig. 12 demonstrates that similar average mole fractions are obtained from NVT-GEMC simulations without and with volume for ∆a ≥ 8.50, and the present simulation data also agree to within the statistical uncertainties with those of Wijmans et al.10 In contrast, when NA , NB is studied in NVT-GEMC simulations without volume moves, then large deviations from the correct results are obtained because due to the lever rule the phases in equilibrium should have different volumes. Introducing volume moves for the nonequimolar system allows one to obtain data in agreement with the simulations for equimolar mixtures. Most importantly, the simulation data for ∆a = 7.80 point to a problem for NVT-GEMC simulations without volume moves even for equimolar systems. Although the results from the present work using this simulation protocol agree with those of Wijmans et al. (see Fig. 12), the NVT-GEMC with volume moves and N PT-GEMC simulations indicate that this ∆a value falls outside of the two-phase region (i.e., below the critical ∆a or above the critical temperature). The mole fraction distributions (see Fig. 13) indicate distinct differences for NVT-GEMC simulations of equimolar systems without and with volumes already for ∆a = 10.0, i.e., far away from the critical point. Introducing volume moves results in broader and more asymmetric mole fraction distributions for ∆a ≤ 10.0, with the differences becoming larger as ∆a is reduced. For the system with ∆a = 7.80, the NVT-GEMC simulation without volume moves still exhibits a clear phase split with bimodal mole fraction distribution (but non-zero probabilities for x A = 0.5), whereas the NVT-GEMC simulation with volume yields a single peak albeit with broad shoulders (see Fig. 13). Of course, it must be noted that the sensitivity to inclusion of volume moves for equimolar systems is a finite-size effect.

J. Chem. Phys. 142, 044902 (2015)

FIG. 13. Liquid–liquid composition histograms computed using different GEMC protocols for a mixture consisting of two different DPD particles with equal size. Data are shown only for systems with NA = NB. The black, red, blue, and green colors indicate the values of the additional repulsive cross-interaction parameter (∆a = 12.5, 10.0, 8.5, and 7.8, respectively) for a system with NA = 500 particles, and the magenta and cyan colors indicate system sizes with NA = 2000 and 4000 particles, where the peak positions are similar as those for the smaller system with a given ∆a. The solid, dashed, and dotted lines denote simulations in the N V T ensemble without and with volume moves, and in the N PT ensemble (only for NA = 500).

To investigate this issue further, additional simulations with a total of 4000 and 8000 particles were carried out. For these larger systems, the simulations at ∆a = 7.80 yield only a single peak even for the protocol without volume moves. However, the distributions for the 4000-particle system remain quite different with the simulation without volume moves yielding a broad peak with a relatively flat top, i.e., a composition split is observed for a significant fraction of the configurations, whereas a much narrower peak is found for the simulations with volume moves. In contrast, the distributions for protocols without and with volumes moves are quite similar for the 8000-particle system at ∆ = 7.80. When probing ∆ = 8.5 in the two-phase region, then one observes the usual finite-size dependence of the composition fluctuations with sharper peaks for larger systems. Nevertheless, even for these larger systems, there remains a minor difference between the simulation protocols where the simulations without volume moves yield sharper peaks (obviously because part of the fluctuations is suppressed), but there are no statistically significant differences in the average values. Comparing the mole fraction distributions from NVTGEMC simulations with volume moves to those from N PTGEMC simulations, one observes slight differences with the distribution for the latter shifted very slightly to a larger miscibility gap at ∆a = 10.0 and 8.5. These differences are caused by slight pressure variations resulting from the larger repulsive cross-interaction parameter that yields a positive excess volume of mixing. In summary, these results indicate that volume moves should be included for NVT-GEMC simulations even for equimolar systems consisting of equal-sized particles and that one should choose the ensemble that best matches the conditions of the actual experimental system for which one wants to explore via DPD simulations.


044902-11


G. DPD simulations for acetic anhydride/(hexane or octane) mixtures

We have thus far developed the rGW-EOS for DPD fluids, generalized the relationship between the Flory-Huggins χ parameter (more specifically, the difference in excess chemical potentials) and the repulsive cross-interaction parameter aij to DPD mixtures with particles of different sizes, and established correct GEMC simulation protocols applicable to LLE calculations for DPD-particle mixtures. To assess the performance of these developments for actual systems, we calculated the LLE for two binary mixtures and compared them to experiment. We chose LLE computations for this assessment as these are one of the most stringent tests for predictive modeling. Specifically, we investigated acetic anhydride/n-hexane and acetic anhydride/n-octane mixtures , where each molecule is represented by a single-bead DPD particle. The intermolecular interactions are similar for these two mixtures, but comparing alkanes of different molecular volumes allowed us to probe the DPD-bead size ratios that are accessible with the new developments. The COSMOtherm software47,48 was used to estimate the infinite-dilution transfer free energies for transferring the alkane molecules into acetic anhydride and vice versa over a range of temperatures and to estimate the molecular volume at T = 298 K that was assumed to be independent of temperature (see supplementary material54 for details). These quantities are input parameters needed to compute χ from Eq. (18). For the two binary mixtures considered here, the LLE are well described by COSMOtherm (see Fig. 14). It should be noted that alternative approaches, such as simulations with molecular

FIG. 14. Liquid–liquid coexistence curves for acetic anhydride/n-hexane (top) and acetic anhydride/n-octane (bottom). The experimental data49 and COSMOtherm predictions are shown as the black lines and filled diamonds, respectively. Circles and squares denote simulation data obtained from N V T and N PT -GEMC simulations, respectively. The cyan, magenta, and green colors indicate DPD predictions using the molecular volumes of acetic anhydride and of the alkane, and their geometric mean, respectively, for vref .

J. Chem. Phys. 142, 044902 (2015)

mechanics force fields,28 quantitative structure-property relationships, or group contribution based methods, could also be used for the prediction of transfer free energies and molecular volumes. The molecular volume ratios are 1.40 (=220/157) for n-hexane versus acetic anhydride and 1.75 (=275/157) for n-octane versus acetic anhydride (see Table I). To explore the robustness of our approach, three different choices of vref (in Eq. (18)) were used for each mixture: (i) the molecular volume of acetic anhydride, (ii) the molecular volume of the alkane, and (iii) the geometric mean of acetic anhydride and alkane molecular volumes. As usual, we set ρref = 5.00 and aref = 15.0 (for PEOS = 41.9), and the remaining pure-phase number densities and repulsive self-interaction parameters were calculated using Eq. (6) and the rGW-EOS, respectively. As described above, we set aij to specific values and find the corresponding χ values and temperatures from Eqs. (19) and (18), respectively. The resulting simulation parameters are summarized in Table I. It should be noted that only for the acetic anhydride/n-hexane mixture with vref = vACE and vref = (vace vhex)1/2 do we obtain pure-phase number densities for both molecules that fall within the 3.00 ≤ ρ0 ≤ 6.00 range considered in the development of Eqs. (20) and (22). In the other cases, either ρ0ace > 6.00 or ρ0oct < 3.00. In contrast, the repulsive cross-interaction parameter always falls within the lower third of the range considered for the parameterization. The LLE obtained from NVT and N PT GEMC simulations with the DPD parameters from Table I are shown in Fig. 14. Although the DPD simulations yield liquid–liquid coexistence curves that are shifted to somewhat lower temperatures than the experimental data,49 an underprediction of the upper critical solution temperature by only 10% is quite remarkable. One should note that this shift in temperature corresponds to a change in the repulsive cross interaction parameter by only one DPD unit. The second remarkable feature is that simulations are in very close agreement for the three vref values considered for each mixture, i.e., even when some of the ρ0 values are outside the parameterization range, Eqs. (20) and (22) appear to hold reasonably well. In addition, this close agreement does not seem to deteriorate when considering the acetic anhydride mixture with the larger volume ratio. As the upper critical solution temperature is approached, the difference between predictions from NVT and N PT simulations increases to a noticeable extent. Different from the athermal mixtures, we observed a significant pressure increase in the NVT simulations for the mutually solvated systems with Psim ranging from about 42.5 at the lowest temperature to about 44.5 at the highest temperature (with PEOS = 41.9). This pressure increase is somewhat smaller for vref = vace and somewhat larger for vref = valk; the reason for this is that the rGW-EOS holds slightly better at higher ρ. In general, if χ > 0 for a NVT simulation, then the pressure of a mixture is expected to be greater than that of the underlying singlecomponent fluids or that of an athermal mixture due to the additional repulsive cross interactions between the beads. With increasing temperature, the mutual miscibilities increase and two competing effects influence the pressure. First, χ and consequently, the repulsive cross interactions decrease; an effect that should mitigate an increase in the pressure. Second,


044902-12


J. Chem. Phys. 142, 044902 (2015)

TABLE I. Simulation parameters for acetic anhydride/n-hexane (A/H) and acetic anhydride/n-octane (A/O) mixtures (ρ ref = 5.00, a ref = 15.0, and PEOS = 41.9). The cross-interaction parameters are provided for the athermal mixture and for the lowest and highest temperatures considered for each system. vref (Å3)

0 ρA

0 ρH

a AA

a HH

a AH(χ = 0)

a AH(Tlow)

a AH(Thigh)

157 220 188

5.00 6.98 5.99

3.58 5.00 4.29

15.0 7.13 10.0

31.6 15.0 21.1

21.9 10.4 15.0

28.0 (224 K) 14.5 (210 K) 19.9 (222 K)

26.4 (279 K) 13.1 (287 K) 18.5 (285 K)

vref (Å3)

0 ρA

0 ρO

a AA

a OO

a AO(χ = 0)

a AO(Tlow)

a AO(Thigh)

157 275 208

5.00 8.74 6.61

2.86 5.00 3.78

15.0 4.28 8.05

53.0 15.0 28.0

28.7 8.03 15.1

35.3 (247 K) 11.7 (227 K) 20.0 (238 K)

33.9 (286 K) 10.6 (293 K) 18.6 (300 K)

the fraction of unfavorable contacts between the two molecular types increases with temperature. The latter appears to be more important for the two systems investigated here because an overall pressure increase is observed. The larger differences in the pressure (by ≈5%) at the higher temperature, that are of course absent in N PT simulations, are likely responsible for the larger differences in the coexistence compositions for NVT and N PT simulations. To eliminate these pressure effects (that would have been much larger for NVT-GEMC simulations for mixtures of differently sized beads with self interaction parameters obtained using the GW-EOS instead of the rGWEOS), we recommend using constant-pressure protocols for DPD simulations of thermal mixtures. Using the pressure as thermodynamic constraint is also more consistent with experiments that are typically conducted at constant pressure. Another observation is that the DPD simulations exhibit a larger miscibility of acetic anhydride in the alkane-rich phases than vice versa, whereas the COSMOtherm predictions show the opposite behavior. At present, it is not clear what may cause these differences. The available experimental data49 show a nearly symmetric LLE curve for the acetic anhydride/n-hexane mixture, whereas that for the acetic anhydride/n-octane is slightly skewed to higher acetic anhydride mole fractions as also predicted by the DPD simulations.

equimolar, symmetric systems, but we recommend to use constant-pressure simulations as much as possible. When applying the new DPD formalism (i.e., the rGW-EOS and the generalized χ-aij relationship) to two actual mixtures, we find very satisfactory predictions of the liquid–liquid coexistence curves for different choices of the reference volume that hold even when the molecular volumes differ by as much as a factor of 1.75. In future work, we plan to explore extensions of this approach toward modeling systems with larger size differences that will require multi-bead representations at least for the larger molecules. ACKNOWLEDGMENTS

Financial support from the National Science Foundation (CBET-1159837) is gratefully acknowledged. Computer resources were provided by the Minnesota Supercomputing Institute. The authors would like to thank Michael Schmidt (P&G) for running the COSMOtherm calculations used in analyzing the liquid–liquid equilibria. P.H.K. acknowledges helpful discussion with Pierre Verstraete (P&G) on using COSMOtherm as a basis for DPD parameterization. 1P. J. Hoogerbrugge and J. M. V. A. Koelman, Europhys. Lett. 19, 155 (1992). 2J. M. V. A. Koelman and P. J. Hoogerbrugge, Europhys. Lett. 21, 363 (1993). 3P.

Español and P. B. Warren, Europhys. Lett. 30, 191 (1995). D. Groot and P. B. Warren, J. Chem. Phys. 107, 4423 (1997). 5S. M. Willemsen, T. J. H. Vlugt, H. C. J. Hoefsloot, and B. Smit, J. Comput. Phys. 147, 507 (1998). 6R. D. Groot and T. J. Madden, J. Chem. Phys. 108, 8713 (1998). 7R. D. Groot, T. J. Madden, and D. J. Tildesley, J. Chem. Phys. 110, 9739 (1999). 8R. D. Groot, Langmuir 16, 7493 (2000). 9R. D. Groot and K. L. Rabone, Biophys. J. 81, 725 (2001). 10C. M. Wijmans, B. Smit, and R. D. Groot, J. Chem. Phys. 114, 7644 (2001). 11S. Yamamoto, Y. Maruyama, and S. Hyodo, J. Chem. Phys. 116, 5842 (2002). 12J. C. Shillcock and R. Lipowsky, J. Chem. Phys. 117, 5048 (2002). 13G. Pan and C. W. Manke, J. Rheol. 46, 1221 (2002). 14L. Rekvig, M. Kranenburg, J. Vreede, B. Hafskjold, and B. Smit, Langmuir 19, 8195 (2003). 15X. Guerrault, B. Rousseau, and J. Farago, J. Chem. Phys. 121, 6538 (2004). 16M. Kranenburg and B. Smit, J. Phys. Chem. B 109, 6553 (2005). 17X. Cao, G. Xu, Y. Li, and Z. Zhang, J. Phys. Chem. A 109, 10418 (2005). 18V. Ortiz, S. O. Nielsen, D. E. Discher, M. L. Klein, R. Lipowsky, and J. Shillcock, J. Phys. Chem. B 109, 17708 (2005). 19J. C. Shillcock and R. Lipowsky, Nat. Mater. 4, 225 (2005). 20J. A. Backer, C. P. Lowe, H. C. J. Hoefsloot, and P. D. Iedema, J. Chem. Phys. 123, 114905 (2005). 21I. V. Pivkin and G. E. Karniadakis, J. Chem. Phys. 124, 184101 (2006). 4R.

V. CONCLUSION

DPD simulations with soft potentials are much better suited than conventional all-atom or coarse-grained MD methods when studying soft-matter at large length and time scales.50–53 In this work, we present three important developments that enable to improve the fidelity of DPD simulations for “thermal” mixtures consisting of chemically dissimilar particles with different pure-phase particle number densities. The first is the development of the rGW-EOS for singlebead DPD fluids (Eq. (14)) that yields pressures with mean unsigned percentage errors that are less than 1% for ρ > 3.00 and an order of magnitude smaller than for the GW-EOS.4 Second, we generalized the relationship between the FloryHuggins χ parameter and the DPD cross-interaction parameter to mixtures with particles possessing different pure-phase number densities (see Eqs. (18)–(22)). Third, we investigated GEMC protocols for DPD simulation and found that volume moves should also be included even for simulations of


044902-13 22V.


Symeonidis, G. E. Karniadakis, and B. Caswell, J. Chem. Phys. 125, 184902 (2006). 23L. Gao, J. Shillcock, and R. Lipowsky, J. Chem. Phys. 126, 015101 (2007). 24J. R. Spaeth, T. Dale, I. G. Kevrekidis, and A. Z. Panagiotopoulos, Ind. Eng. Chem. Res. 50, 69 (2011). 25G. Kacar, E. A. J. F. Peters, and G. de With, Europhys. Lett. 102, 40009 (2013). 26A. Vishnyakov, M. Lee, and A. V. Neimark, J. Phys. Chem. Lett. 4, 797 (2013). 27M. Lee, A. Vishnyakov, and A. V. Neimark, J. Phys. Chem. B 117, 10304 (2013). 28E. Mayoral and A. G. Goicochea, J. Chem. Phys. 138, 094703 (2013). 29I. Pagonabarraga and D. Frenkel, J. Chem. Phys. 115, 5015 (2001). 30P. B. Warren, Phys. Rev. E 68, 066702 (2003). 31S. W. Smith, C. K. Hall, and B. D. Freeman, J. Comput. Phys. 134, 16 (1997). 32A. J. Schultz, C. K. Hall, and J. Genzer, Macromolecules 38, 3007 (2005). 33J. Liu, T. L. Bowman, and J. R. Elliott, Ind. Eng. Chem. Res. 33, 957 (1994). 34D. A. McQuarrie, Statistical Mechanics (Harper Collins, New York, 1976). 35A. Maiti and S. McGrother, J. Chem. Phys. 120, 1594 (2004). 36K. P. Travis, M. Bankhead, K. Good, and S. L. Owens, J. Chem. Phys. 127, 014109 (2007). 37See http://www.chem.umn.edu/groups/siepmann/software.html. for Monte Carlo for complex chemical systems-Minnesota. 38N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).

J. Chem. Phys. 142, 044902 (2015) 39W.

W. Wood and F. R. Parker, J. Chem. Phys. 27, 720 (1957). R. McDonald, Mol. Phys. 23, 41 (1972). 41A. Panagiotopoulos, N. Quirke, M. Stapleton, and D. J. Tildesley, Mol. Phys. 63, 527 (1988). 42J. I. Siepmann and D. Frenkel, Mol. Phys. 75, 59 (1992). 43K. Esselink, L. D. J. C. Loyens, and B. Smit, Phys. Rev. E 51, 1560 (1995). 44M. G. Martin and J. I. Siepmann, Theor. Chem. Acc. 99, 347 (1998). 45A. Ben-Naim, Statistical Thermodynamics for Chemists and Biochemists (Plenum, New York, 1992). 46A. Z. Panagiotopoulos, Mol. Phys. 61, 813 (1987). 47A. Klamt, V. Jonas, T. Bürger, and J. C. W. Lohrenz, J. Phys. Chem. A 102, 5074 (1998). 48A. Klamt and F. Eckert, Fluid Phase Equilib. 172, 43 (2000). 49M. Aboy, S. Villa, N. Riesco, J. A. González, I. G. Fuente, and J. C. Cobos, J. Chem. Eng. Data 47, 950 (2002). 50S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. de Vries, J. Phys. Chem. B 111, 7812 (2007). 51L. Monticelli, S. K. Kandasamy, X. Periole, R. G. Larson, D. P. Tieleman, and S. J. Marrink, J. Chem. Theory Comput. 4, 819 (2008). 52S. J. Marrink and D. P. Tieleman, Chem. Soc. Rev. 42, 6801 (2013). 53W. Shinoda, R. Devane, and M. L. Klein, Mol. Simul. 33, 27 (2007). 54See supplementary material at http://dx.doi.org/10.1063/1.4905918 for tables listing the numerical data for pressures of pure phases and mixtures, radial distribution functions, infinite dilution excess chemical potentials, and coexistence mole fraction data for LLE simulations. 40I.


A polarizable coarse-grained water model for dissipative particle dynamics.

Dissipative particle dynamics with an effective pair potential from integral equation theory of molecular liquids.

Parametric study of particle sedimentation by dissipative particle dynamics simulation.

Distortion and flow of nematics simulated by dissipative particle dynamics.

Equation of state and adsorption dynamics of soft microgel particles at an air-water interface.

Multiscale simulation of ideal mixtures using smoothed dissipative particle dynamics.

Multiscale Modeling of Flow Induced Thrombogenicity With Dissipative Particle Dynamics and Molecular Dynamics.

Efficient Schmidt number scaling in dissipative particle dynamics.

Co-micellization behavior in poloxamers: dissipative particle dynamics study.

Modeling nanoscale hydrodynamics by smoothed dissipative particle dynamics.

Hybrid molecular-continuum simulations using smoothed dissipative particle dynamics.

Emergence of glasslike dynamics for dissipative and strongly interacting bosons.

An efficient dissipative particle dynamics-based algorithm for simulating electrolyte solutions.

Generalized mapping of multi-body dissipative particle dynamics onto fluid compressibility and the Flory-Huggins theory.

Effects of Varying Particle Sizes and Different Types of LDH-Modified Anthracite in Simulated Test Columns for Phosphorous Removal.

SIMS instrumentation and methodology for mapping of co-localized molecules.

Equation of State for Phospholipid Self-Assembly.

In vitro and in vivo evaluation of SN-38 nanocrystals with different particle sizes.

Carbon mineralization in two ultisols amended with different sources and particle sizes of pyrolyzed biochar.

Determination of macroscopic transport coefficients of a dissipative particle dynamics solvent.

Preparation and evaluation of tilmicosin-loaded hydrogenated castor oil nanoparticle suspensions of different particle sizes.

Dissipative particle dynamics simulation study of poly(2-oxazoline)-based multicompartment micelle nanoreactor.

Genotoxicity assessment of magnetic iron oxide nanoparticles with different particle sizes and surface coatings.

Dissipative particle dynamics of triblock copolymer melts: a midblock conformational study at moderate segregation.