A chain-revised Groot-Warren equation of state (crGW-EOS) was developed and tested to describe systems of homo-oligomeric chains in the framework of dissipative particle dynamics (DPD). First, thermodynamic perturbation theory is applied to introduce correction terms that account for the reduction in pressure with an increasing number of bonds at constant bead number density. Then, this EOS is modified by introducing a set of switching functions that yields an accurate second virial coefficient in the low-density limit. The crGW-EOS offers several improvements over the revised Groot-Warren equation of state and Groot-Warren equation of state for chain molecules. We tested the crGW-EOS by using it to predict the pressure of oligomeric systems and the B2 virial coefficient of chain DPD particles for a range of bond lengths. Additionally, a method is developed for determining the strength of cross-interaction parameters between chains of different compositions and sizes and for thermal and athermal mixtures. We explored how different levels of coarse-graining affect the upper-critical solution temperature.

In 1993, Hoogerbrugge and Koelman developed the dissipative particle dynamics (DPD) method.1,2 DPD is a coarse-graining method used to model mesoscopic phenomena of complex fluids by aggregating several atoms or small molecules together as a single particle whose interaction with other particles is described by soft, purely repulsive conservative potentials. Unlike other coarse-graining approaches such as the MARTINI,3 Shinoda-DeVane-Klein,4 statistical associating fluid theory-γ Mie,5 and transferable potentials for phase equilibria–coarse-grain6 force fields, DPD enhances simulation speed by using a softer potential (allowing for a longer time step in molecular dynamics) and neglecting long-range interactions (reducing the number of distance calculations and/or force evaluations). In the DPD method, particles are assigned random and dissipative forces that can be tuned to satisfy the fluctuation dissipation theorem and to maintain a constant temperature in the simulation.7,8 These forces can also be tuned with the time step size to yield correct hydrodynamical behavior, ensure conservation of momentum, and follow the equipartition theorem.8 

In 1997, Groot and Warren made the DPD method more useful by introducing the Groot-Warren equation of state (GW-EOS).7 Groot and Warren developed the GW-EOS by using the virial theorem and fitting simulation data to relate the force field repulsive constant, aij, to the thermodynamic compressibility in systems comprising particles that interact with the DPD potential.7 The GW-EOS is given by

pGW=ρi0kBT+αaiiρi02,
(1)

where pGW is the pressure of the system, T is the absolute temperature, kB is the Boltzmann constant, aii is the repulsive force constant, and ρi is the pure-phase number density of particle type i in the fluid where the pure-phase density is defined in relation to a reference system by a method developed by Kacar et al.9 Groot and Warren found α = 0.101 ± 0.001 through fitting to systems with three different values of aii over a range of densities.10 With the GW-EOS, a new method to connect interaction strength between beads of the same type to the macroscopic properties of a fluid was introduced.7 The GW-EOS framework has allowed researchers to explore liquid-liquid equilibria (LLE), interfacial tensions, the formation of surfactant mesophases, homo- and block polymer melts, and many other self-assembly and phase-separation phenomena.11–29 

But the GW-EOS method has certain drawbacks. First, the GW-EOS relies on the Flory-Huggins theory to parametrize interaction strengths between different bead types in mixtures.9,10,30 Flory-Huggins theory assumes that all particles are of the same size; in the DPD framework, this means that each particle type is required to have the same pure-phase density (i.e., the same size and self-interaction strength) regardless of the actual size and chemistry of the fragments being represented in DPD.30 Therefore, the GW-EOS cannot be used to accurately predict the thermodynamic and mechanical properties of mixtures.9,10,30 Another drawback is that the GW-EOS framework yields inaccurate predictions of the pressure of a system as density is decreased (i.e., pressure is decreased). Groot and Warren were aware of this, noting that inaccuracy occurs because the leading coefficient in the quadratic density term of the pressure equation of state is approximated as a constant for systems with a sufficiently high density.7 The result is that the GW-EOS improperly recovers virial coefficients in the low-density limit.7 Maiti and McGrother noted that the value of α is not constant, but rather a weak function of aii and ρ.10 Similarly, Liyana-Arachchi et al. demonstrated that Groot and Warren’s approximation of α deviates significantly as density decreases.31 

Despite GW-EOS’s shortcomings, the popularity of the method has caused researchers to refine it. For instance, by building on the work of Maiti and McGrother,10 Travis et al. united the relation for the free energy of a DPD fluid with Scatchard-Hildebrand regular solution theory, allowing the interaction between particles to be different in both size and chemistry.30 With their approach, the self-interaction strengths between like beads were parametrized with solubility parameters, and cross-interactions were parametrized with differences in cohesive energy density between each bead type.30 

To create a more applicable and accurate equation of state, particularly at lower particle densities (allowing for a smaller number of DPD particles to represent a given system), Liyana-Arachchi et al. developed the revised Groot-Warren equation of state (rGW-EOS).31 The rGW-EOS allows researchers to introduce the B2 virial coefficient and a set of switching functions to restore virial characteristics expected at low densities31 and to adjust how much the quadratic terms contribute to the pressure in the equation of state. The rGW-EOS is given by31 

prGW=ρi0kBT+[B2(aii)FB2(ρi0)+aiiFα(ρi0)](ρi0)2,
(2)

with

FB2ρ=11+ρ3
(3)

and

Fαρ=c1ρ21+c2ρ2,
(4)

where B2(aii) is the second virial coefficient of a unary DPD fluid with the repulsive force constant of aii, c1 and c2 are equal to 0.0802 and 0.7787, respectively (determined by fitting with test simulation data). The rGW-EOS expands the range of values over which the equation of state is accurate for single bead representations, and it improves the overall accuracy of the equation of state approach.31 When using the rGW-EOS to analyze a range of state points similar to those analyzed by Groot and Warren,7,31 the mean of unsigned percent error falls to 0.7% from the 11.5% yielded by GW-EOS. Even over a range of state points common to DPD simulations, the rGW-EOS is far more accurate than the GW-EOS.

When the rGW-EOS and GW-EOS are applied to systems composed of chain molecules, it is relevant whether one should choose the bead density of the system or the chain density of the system as the input density. For systems comprising single-bead chains, there is no difference between the bead and chain densities of the system. However, for a system comprising multi-bead chains, the difference between the bead and chain densities becomes significant. Typical DPD systems are often packed tightly enough that researchers make an assumption that the bead density can be used as the input density. At low densities, this leads to poor parametrization as the ideal gas and virial terms (both largely dependent on the chain density) begin to have a dominant contribution to the pressure in the equation of state. To circumvent this issue, one can use the bead density while making the equation of state sensitive to the number of beads in each chain. However, because the rGW-EOS and GW-EOS models are insensitive to the bonding details of the chains—such as the bond length—this approach still leads to the incorrect calculation of appropriate parameters for each system.

Thus, this work is an extension of the rGW-EOS in the work of Liyana-Arachchi et al. To address the need for the rGW-EOS to describe homo-oligomer chains, we devise a new equation of state called the chain-revised Groot-Warren equation of state (crGW-EOS). It is developed using first-order thermodynamic perturbation theory (TPT1) and Wertheim’s theory of polymerization, which together account for the reduction in pressure when bonds are introduced.32–35 We also validate the crGW-EOS through a series of Monte Carlo (MC) simulations. Finally, we include the effects of bonding in thermal and athermal mixtures on the rGW-EOS framework.31 

Within the framework of DPD, the interactions between particles are described by a soft and purely repulsive conservative potential

urij=aij21rijrcut2rijrcut10rijrcut>1,
(5)

where aij describes the strength of the interaction between two particles, i and j, which are separated by a distance, rij. This potential is truncated at a distance in dimensionless DPD units, rcut, where it falls to 0. In this work, a harmonic potential between bonded beads (1-2 interactions) is used to model a number, m, of DPD particles chained together in homo-oligomer chains. This potential has a coarse-grained spring constant, kspring, and a corresponding bond length, l = rbond/rcut. The strength of the bonds does not affect the system pressure if the bonds have a rigidity kspring > 100. At bond lengths l < 0.25, intra-chain (1-3) interactions cause issues in properly sampling the bond length distribution; because of this, a larger spring constant is needed to properly constrain the average simulated bond length.

All simulations in this work are done using the MCCCS-MN (Monte Carlo for Complex Chemical Systems–Minnesota) software suite. The kd-tree data structure for computing particle-particle interactions is used to add significant efficiency.36 Unless otherwise noted, probabilities for MC moves (e.g., center-of-mass translation, center-of-mass rotation, and configurational bias) in each simulation are chosen to be equal to the relative number of degrees of freedom involved in each move type.

To tune fitting parameters η and β (note: these parameters are explained in detail in Sec. III A 1) for the crGW-EOS and to validate the crGW-EOS against the GW-EOS and rGW-EOS, we perform several MC simulations in a one-box NVT ensemble using a wide range of system parameters. These include the number of DPD beads bonded in each chain, m (2, 4, or 8), bead density, ρ (3, 4, 5, 6, or 7), bond length, l (0.3, 0.4, 0.5, or 0.6), and bead-bead self-interaction parameter, aii (10, 20, 30, or 40). In each simulation, 3072 DPD beads are modeled, which corresponds to 1536, 738, and 384 molecules when modeled in 2-, 4-, or 8-bead chains, respectively. The volume of the simulation box is chosen to simulate the correct bead density for each simulation. The results are averaged across four independent simulations.

The second virial coefficient for chain molecules, B2chain, is calculated based on a method used by Martin and Siepmann.37 Martin and Siepmann’s method involves calculating the second virial coefficient of a fully flexible molecule,38,39

B2(T)=2π0eUinter(r12)kBTα1,α21r122dr12,
(6)

where Uinter is the potential between two molecules whose centers of mass are a distance, r12, apart. The α1,α2 represents the canonical ensemble averaged over configurations weighted solely on their Boltzmann weights. The configurations are sampled in simulations comprising two boxes, each box containing a single chain. The conformational space is sampled through rotational and configurational moves. Center of mass separations larger than twice the maximum extent of a chain are not considered. The intermolecular energy is calculated for r12 values in step sizes of 0.05 Å/rcut. We run these two-box simulations using the same range of system parameters (m, l, and aii) as with the one-box tuning simulations. The results are averaged across eight independent simulations.

To validate a cross-interaction term, aij, added to the crGW-EOS, we run MC simulations of mixtures in the NpT and NVT ensembles. First, three fictitious mixtures are tested in the athermal regime: a 1-bead fluid, Vh1, with a molecular volume similar to that of n-hexane (220 Å3) in a mixture with a 2-bead fluid, Vh2, of the same volume (220 Å3); Vh1 in a mixture with a 2-bead fluid, Vo2, with a molecular volume similar to that of n-octane (275 Å3); and Vh2 in a mixture with Vo2 (these molecular volumes were estimated by Liyana-Arachchi et al.31 using the COSMOTherm software).40,41 For parametrization purposes, we use a reference volume similar to the geometric mean between the volume of Vh1 and the volume of a single bead in Vh2, 150 Å3. This density also corresponds to five water molecules per reference bead. The water reference fluid has a density, ρref, of 5.0, which yields a cutoff distance, rcut, of 9.086 and a reference pressure, pref, of 41.89. The pure-phase density, ρ0, of the fictitious n-hexane and n-octane fluids is set equal to 3.432 and 2.780, respectively, to match experimental conditions. The bond lengths, l, for Vh2 and Vo2 are set to 0.39 and 0.5, respectively, to match distances between centers of mass obtained from united atom simulations. Table I shows the self- and cross-interaction parameters for each mixture and each equation of state. Each simulation includes 2000 molecules, and each component’s mole fraction ranges from xi = 0.0 to xi = 1.0 in step sizes of 0.1. The simulations have an equilibration period of at least 70 000 MC cycles followed by a production period of at least 70 000 MC cycles averaged across eight independent trajectories.

TABLE I.

The self- and cross-interaction parameters for each athermal mixture using the crGW-EOS, rGW-EOS, or GW-EOS. For the rGW-EOS and the GW-EOS, the bead density was used as an input.

aii
ModelVh1Vh2aij
crGW-EOS 35.0 8.57 17.3 
rGW-EOS 35.0 7.43 16.1 
GW-EOS 32.3 7.36 15.4 
 Vh1 Vo2  
crGW-EOS 35.0 13.2 21.5 
rGW-EOS 35.0 11.9 20.4 
GW-EOS 32.3 11.6 19.4 
 Vh2 Vo2  
crGW-EOS 8.57 13.2 10.6 
rGW-EOS 7.46 11.9 9.4 
GW-EOS 7.36 11.6 9.24 
aii
ModelVh1Vh2aij
crGW-EOS 35.0 8.57 17.3 
rGW-EOS 35.0 7.43 16.1 
GW-EOS 32.3 7.36 15.4 
 Vh1 Vo2  
crGW-EOS 35.0 13.2 21.5 
rGW-EOS 35.0 11.9 20.4 
GW-EOS 32.3 11.6 19.4 
 Vh2 Vo2  
crGW-EOS 8.57 13.2 10.6 
rGW-EOS 7.46 11.9 9.4 
GW-EOS 7.36 11.6 9.24 

To calculate a Flory-Huggins χ-parameter, we first calculate 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,42,43

ΔGi,jtrans=μi,jex=RTlnρliq,jiρvapi,
(7)

where ρliq,ji and ρvapi are the average number densities of particle i in the liquid and vapor phases, respectively, and R is the molar gas constant. We follow the work of Liyana-Arachchi et al. in performing MC simulations in a two-box pseudo-osmotic Gibbs ensemble. In these simulations, volume moves are only performed on the liquid box. The solute, i, is the only particle allowed to swap between the boxes. The liquid phase is kept at a constant pressure of 41.89, based on a reference water bead which comprises three water molecules and which has a reference density of 5.0. In this work, the solvent is made of monomer units, and we test four different levels of coarse-graining for the solute particle: a dimer, trimer, tetramer, and hexamer (where each “monomer” unit is represented by a bead). The number of solvent chains in each simulation (500, 500, 970, and 3000 for the dimer, trimer, tetramer, and hexamer simulations, respectively) differs to ensure that the liquid box is large enough to accommodate the larger solute molecules. The bond length of each chain ranges between 0.1 and 1.5 in step sizes of 0.1 (an additional bond length of 2.5 was also tested). The coarse-grained pure-phase solvent densities range from 2.0 to 5.5, 2.0 to 4.5, 2.0 to 4.5, and 2.0 to 3.5 for the dimer, trimer, tetramer, and hexamer systems, respectively, in step sizes of 0.5. These are shown in Table II, along with the corresponding density simulated and the solvent-solvent interaction parameter calculated at each density. Note that the solute-solute interaction parameter is set to 0 because there is only one solute in each simulation. The cross-interaction parameters, chosen to be similar to athermal mixing values, are shown in Table III. Volume moves are only allowed on the liquid box, and a biasing potential is introduced with swap moves to ensure that the solute is sampled evenly in both boxes. The equilibration is run for at least 30 000 MC cycles, and the production is also run for at least 30 000 MC cycles with four independent trajectories.

TABLE II.

The pure-phase density of each solvent fluid. The rGW-EOS is used to generate the ajj parameters. The left-most column is the input pure-phase density for each system, and the middle column is the simulated density of a neat solvent system generated from NpT Gibbs ensemble simulations to show the deviation from the value predicted by the equation of state.

Solvent parameters
Input densitySimulated densityajj
2.0 1.961 126.6 
2.5 2.442 73.25 
3.0 2.962 47.72 
3.5 3.4813 33.52 
4.0 3.994 24.80 
4.5 4.505 19.05 
5.0 5.016 15.07 
5.5 5.516 12.19 
Solvent parameters
Input densitySimulated densityajj
2.0 1.961 126.6 
2.5 2.442 73.25 
3.0 2.962 47.72 
3.5 3.4813 33.52 
4.0 3.994 24.80 
4.5 4.505 19.05 
5.0 5.016 15.07 
5.5 5.516 12.19 
TABLE III.

: The cross-interaction parameters between the solute and solvent beads. For each solute, aij is chosen to be similar to an athermal mixing value.

Solute typeaij
Dimer 23.39 
Trimer 15.74 
Tetramer 17.73 
Hexamer 7.33 
Solute typeaij
Dimer 23.39 
Trimer 15.74 
Tetramer 17.73 
Hexamer 7.33 

To evaluate which equation of state—crGW-EOS, rGW-EOS, or GW-EOS—gives the most reasonable upper-critical solution temperature (UCST), three fictitious mixtures are tested for each equation of state with MC simulations in a two-box NpT Gibbs ensemble: Vh1 with Vh2, Vh1 with Vo2, and Vh2 with Vo2. In each simulation, there are 1000 chains of each molecule type, and the system is initialized with only one molecule type per box. The χ−1 values we test range from 0.2 to 0.56; this range is selected based on cross-interaction parameters that lead to mixtures being neither totally miscible nor immiscible. The self-interaction parameters are parametrized to correspond to liquid n-hexane or n-octane (i.e., they are calculated from the number densities of liquid n-hexane or n-octane) at 273 K and 1 atm. The self-interaction parameters and bond lengths are outlined in Table IV. The DPD reference fluid is taken as a bead with a density of 5.0 (corresponding to five water molecules per reference bead) and pressure of 41.89, which yields an rcut of 9.09. The systems undergo equilibration for 180 000 total MC cycles (40 000 of which are just volume equilibration and the remaining involve swap moves). After equilibration, four independent simulations are spawned for production, which are run for at least 70 000 MC cycles.

TABLE IV.

Simulation parameters for octane and hexane using each equation of state method.

Modelm-beadsaiili
crGW (hexane) 35.0 … 
 8.57 0.39 
crGW (octane) 13.2 0.5 
rGW (hexane) 35.0 … 
 7.43 0.39 
rGW (octane) 11.9 0.5 
GW (hexane) 32.3 … 
 7.36 0.39 
GW (octane) 11.6 0.5 
Pure-phase density (hexane) 3.432   
Pure-phase density (octane) 2.780   
Modelm-beadsaiili
crGW (hexane) 35.0 … 
 8.57 0.39 
crGW (octane) 13.2 0.5 
rGW (hexane) 35.0 … 
 7.43 0.39 
rGW (octane) 11.9 0.5 
GW (hexane) 32.3 … 
 7.36 0.39 
GW (octane) 11.6 0.5 
Pure-phase density (hexane) 3.432   
Pure-phase density (octane) 2.780   

To test the crGW-EOS on a real system, we study mixtures of n-hexane/acetic anhydride and n-octane/acetic anhydride using MC simulations in a two-box NpT Gibbs ensemble. These mixtures are parameterized with the COSMOTherm software. Prior to these simulations, harmonic bond parameters are estimated using gas phase MC simulations in which the molecules are described by the COMPASS (Condensed-phase Optimized Molecular Potentials for Atomistic Simulation Studies) force field.44 Each mixture simulation is parametrized using a reference liquid comprising three water molecules per bead, where ρref = 5.0, vref = 90 Å3, which yields rcut = 7.66 Å, and aii,ref = 15.0, which yields p = 41.89. From these parameterizations, the resulting self-interaction parameters for hexane beads range from 25.74 at 160 K to 25.97 at 400 K. For octane beads, the self-interaction parameters range from 40.47 at 160 K to 41.23 at 400 K. The acetic anhydride self-interaction parameters are 13.18 and 53.73 for the 2-bead and 1-bead representations, respectively (note that these do not depend on temperature). Two sets of simulations are performed; one in which the system is initialized with only one molecule type per box. In the other, each box is initialized with an equimolar mixture of the two molecules. Each mixture simulation has a total system size of 2000 molecules split equally amongst the two molecule types. The mixtures are allowed at least 300 000 MC cycles of equilibration (40 000 of which are just volume equilibration and the remaining involve swap moves). After equilibration, four independent simulations (for each mixture) are spawned for production and ran for at least 100 000 MC cycles.

The equation of state we propose below for bonded molecules requires the cavity correlation function, y(r) of the monomeric DPD fluid and its density derivative of its logarithm (at the bond length of the chain molecule). y(r) is defined as

y(r)=exp[βU(r)]g(r)
(8)
ln[y(r)]=βU(r)+ln[g(r)]
(9)
ρln[y(r)]=ρln[g(r)]=ρw(r)(kBT=1),
(10)

where w(r) is the potential of mean force, g(r) is the radial distribution function of the monomeric DPD fluid, and U(r) is the conservative part of the DPD interaction potential. Instead of using expensive simulations, g(r) for any particular choice of {a, ρ} was calculated numerically (and on the fly) using the Ornstein-Zernike equation with the hypernetted chain closure (OZ-HNC). For the relatively soft DPD potential, the g(r) obtained from OZ-HNC converges rapidly and is highly accurate (confirmed by comparing to explicit simulations). The derivative of g(r) with respect to density at the bond length of interest was calculated numerically from finite differences, i.e., w(r, {a, ρ + δρ}) − w(r, {a, ρδρ})/2δρ, with δρ = 0.05. Warren’s open-source SunlightDPD hypernetted chain solver45 is an alternative DPD-specific OZ-HNC solver that can be used to pre-calculate and tabulate these derivatives ahead of the chain equation of state calculations.

Before deriving the crGW-EOS, we placed a set of constraints. First, the crGW-EOS needed to be consistent with the rGW-EOS framework, both physically and mathematically, with the ability to reduce to the rGW-EOS in the 1-bead per chain limit and to a modified rGW-EOS in the limit of zero bond length and a bond force constant of zero. Second, in the low-density limit, the crGW-EOS had to depend on the chain density rather than the bead density. This would ensure the crGW-EOS’s adherence to thermodynamics in the ideal gas limit. In the high-density limit, the quadratic term from the GW-EOS and rGW-EOS produced better predictions when the bead density was used instead of the chain density; therefore, the quadratic term in density would depend on the overall bead density. Any corrective terms would depend on the properties of the interaction between the chain’s monomer units. Additionally, for improved low-density predictions, we expected the B2 virial coefficient, which we develop in the second half of this section, to play a role in describing the interactions of these systems.

The approach of Johnson et al.35 for the pressure of tangent hard sphere and Lennard Jonesium chains, based on Wertheim’s theory of polymerization,32,33 is a useful starting point that meets some of these criteria. Using TPT1, Johnson et al. derived a relation for the Helmholtz free energy of a chain fluid.34,35 In TPT1, it is assumed that the structure of a dense chain fluid is dominated by short range repulsive interactions so that the structure of a chain fluid is similar to that of a monomeric fluid at the same monomeric density. Therefore, the monomer fluid is chosen as the reference, and the effect of bonds on the structure of the system is assumed to cause small perturbation to the reference system, and the equation of state they derived is the following:

p=pRm1mρRkBTm1mρR2[lnyRl]ρR,
(11)

where pR is the pressure of the monomeric reference fluid, the second term captures the difference in ideal gas pressure between the monomeric fluid vs. the chain fluid (i.e., ρRkBTvs. (ρR/m)kBT) and the final term captures the effect of intermolecular bonding. We can modify this equation systematically to include both the high and low density limits as follows:

1. High density limit

At high densities, the repulsive interactions dominate over the ideal gas and (chain) second virial coefficient contributions to the pressure. We can plug in either the GW-EOS or the more accurate rGW-EOS for the pressure of the reference fluid and focus on the effect of introducing bonds in a DPD fluid. However, we find that replacing the (m − 1)/m prefactor in the third term with η(m − 1)/(m + β), where η and β are fitting parameters, allows for a much better match with the simulation results for bead densities, ρ ∈ (3, 5) and α ∈ (10–50). That is, the equation of state is

p=pRm1mρRkBTηm1m+βρR2[lnyRl]ρR.
(12)
a. Physical interpretation of the yR(r) term.

yR(r) represents the probability that there is a cavity at r given that there is a cavity at the origin. For hard spheres, the size of the cavity needed to insert the hard sphere is well defined—so, yR(r = l) represents the conditional probability of a “dimer” shaped cavity with a bond length l given that there is already a cavity large enough to insert a hard sphere at the origin. For soft solutes, such a description of a cavity is somewhat qualitative. So, the first order approximation for the probability to create a cavity needed to introduce a chain of m DPD beads is not [y2(l)]m1 [Eq. (4) of Johnson et al.35], but slightly different, and the empirically introduced parameters, β and η, attempt to capture this difference. At this point, we cannot justify them more rigorously.

An intuitive way to understand the origin of the ρR2ρRln[yR(l)] term in Eq. (11) can be as follows:

ρR2ρRlnyR(l)=ρR2ρRwR(l).
(13)

Now the potential of mean force, w(r), is a free energy in the canonical ensemble. The pressure is the volume derivative of the Helmholtz free energy, P = −∂A/∂V. Substituting ρR ∼ 1/VR, we get ΔP=ρR2A/ρR=ρR2wR(l)/ρR for each bond (of length l) that is created, with (m − 1) monomer bonds per chain. There is an additional 1/m factor since the number of molecules in chain fluid compared to the number of molecules in the reference monomer fluid is reduced by a factor of m. Consequently, the drop in pressure due to the formation of bonds (in addition to the ideal gas term) is ΔPm1mρR2wR(l)/ρR.

2. Incorporating the low density limit and the role of B2

In the low density limit, we expect the pressure of the chain fluid to behave as P/kBT=ρchain+B2ρchain2, where B2 is the second virial coefficient between two chains. We therefore employ the same strategy as in the development of the rGW-EOS and introduce two switching functions to interpolate between the two limits. The equation of state, that we term the “chain revised GW-EOS” (crGW-EOS) thus becomes

pkBT=ρchain+FB2ρchainB2chainρchain2+Fαmρchainmρchain2×aiiηm1m+β[lnyRl]mρchain,
(14)

where chain = ρR, and FB2 and Fα are the same switching functions present in the rGW-EOS,31 and B2chain is the DPD chain virial coefficient. We describe the calculation of B2chain in Sec. III A 2 a.

a. Calculation of B2chain.

The inter-chain second virial coefficient is a function of the chain length (m), bond length (l), and the bead-bead repulsion parameter (a). Unlike a monomeric fluid, it is not possible to directly integrate the Mayer function and fit to a suitable expression as was done for the development of the rGW-EOS. As described in Sec. II, B2chain=B2chain(a,m,l) was calculated for a range of relevant parameter values. The results were then fit to a 3 parameter equation that satisfied certain physical constraints—(i) we must recover the monomeric virial coefficient in the limit of a 1-bead chain. (ii) As the bond length decreases to 0, beads in the chain will overlap entirely and form a monomer made of many beads. The interaction strength between beads will scale with the number of beads in each chain. (iii) For zero inter bead repulsion, the virial coefficient must reduce to zero. Mathematically, these constraints are

limm1B2chainaii,m,l=B2monoaii,
(15)
liml0B2chainaii,m,l=B2monom2aii,
(16)
limaii0B2chainaii,m,l=0.
(17)

These restrictions suggest the following approximate form for the chain second virial coefficient relation:

B2chainaii,m,l=g(aii,l,m)B2monof  m,laii.
(18)

To find B2mono, the Mayer function for a single bead DPD fluid was computed as

B2monoaii=2π0rc(exp[u(aii,r)/kBT]1)r2dr.
(19)

Upon plugging in the DPD potential [Eq. (5)] and solving this complex integral, we end up with the following expression for B2mono:

B2monoaii=2π1+aiiπ2Erfaii2aii32+eaii22aii13,
(20)

where Erf is the Gauss error function. In contrast to the parametrized relation for B2mono developed along with the rGW-EOS,31 this function saturates to the proper hard sphere limit of 2π3 for a large aii, and it is valid for any positive valued aii. This is important because most of the aii-dependence in B2chainaii,m,l arises from the B2mono factor. Figure 1 shows how our relation for B2mono compares to the one developed with the rGW-EOS. Both equations produce similar values for B2mono for all aii less than 50. At aii > 50, the two relations diverge.

FIG. 1.

A comparison of Eq. (20) with a similar function parametrized by Liyana-Arachchi et al. The inset graph is a zoomed-in view of the smaller aii values. The red diamonds and black circles denote the resulting B2mono values at each aii using the rGW-EOS and crGW-EOS, respectively.

FIG. 1.

A comparison of Eq. (20) with a similar function parametrized by Liyana-Arachchi et al. The inset graph is a zoomed-in view of the smaller aii values. The red diamonds and black circles denote the resulting B2mono values at each aii using the rGW-EOS and crGW-EOS, respectively.

Close modal

To find a functional form for g(aii, l, m) and f  m,l, we had to make them independent of one another. To do this, we found a functional relation for f  m,l because Eq. (16) implied

liml0g(aii,l,m)=1.
(21)

We fit simulation data of chain virial coefficients of chains with bond lengths nearing 0 to obtain the function f  m,l, approximated as

f  m,l=m2l+13.
(22)

This suggested that the monomeric contribution to the chain virial coefficient would saturate quickly as the number of beads in the chain increased. With this relation, we isolated the function g(aii, l, m) by factoring out the contribution from the monomeric portion of the equation

B2chain, sim.aii,m,lB2monof  m,laii=g(aii,l,m).
(23)

The resultant quantity could be described by the approximation

g(aii,l,m)=A1γm2+γm1lδ+1,
(24)

where γ, A, and δ are fitting parameters (note that this approximation’s functional form was known because of its limits when m approaches 1 and l approaches 0). A least squares fit to simulation data yields γ = 0.926, A = 1.954, and δ = 1.57. These parameters yielded accurate predictions for the chain virial coefficient, with a mean unsigned percent error (MUPE) of 1.6%. The relative errors in the predicted B2chain from our equation compared to simulated values (see Fig. S1 of the supplementary material) show no pattern to improve upon this equation; an exact relation would be very complex and lacks a simple mathematical solution.

With B2chain, the crGW-EOS is fully described as in Eq. (14),

pkBT=ρchain+FB2ρchainB2chainaii,m,lρchain2     +Fαmρchainmρchain2aiiηm1m+βlnyRlmρchain.

To fit and validate the full pressure equation of state for a DPD chain fluid, we performed simulations, described in Sec. II B, over a wide range of parameters common to DPD simulations. We fit the data from these simulations to the crGW-EOS using a least squares fit method. Upon fitting, we determined fitting parameters of η = 41.083 and β = 1.967, which yielded a 2.5% MUPE. The relative error absolute value varied up to 11.5%.

Figure 2 shows the mismatch between the predicted pressure and the simulated pressure for a given set of parameters. The predicted pressures are distributed around 0, with a slight tendency to overpredict the simulated pressures. Like in the GW-EOS and rGW-EOS, as the bead density and the pressure decrease, the expected error in predicted pressure from the equation of state increases. Typically, this does not present a problem if the reference volume is chosen with a moderate pure-phase bead density for the molecule. The clearest remaining pattern is the mismatch between simulation and experiment, which arises from the number of beads in the chain. As the number of beads in the chain increases, the pressure for a system of 2-bead chains moves from slight under-prediction to over-prediction. Higher order corrections to the equation of state are needed to further correct for this over-prediction.

FIG. 2.

A depiction of the percent mismatch between simulated pressure and pressure predicted by the crGW-EOS. Each subplot is color- and symbol-coded to different parameters to explore patterns in the pressure mismatch. The top left graph is coded according to the bond length l. The top right graph is coded by the bead density ρ. The bottom left graph is coded to the bead-bead interaction parameter aii. The bottom right graph is coded according to the number of beads in each chain m.

FIG. 2.

A depiction of the percent mismatch between simulated pressure and pressure predicted by the crGW-EOS. Each subplot is color- and symbol-coded to different parameters to explore patterns in the pressure mismatch. The top left graph is coded according to the bond length l. The top right graph is coded by the bead density ρ. The bottom left graph is coded to the bead-bead interaction parameter aii. The bottom right graph is coded according to the number of beads in each chain m.

Close modal

Compared to the GW-EOS and rGW-EOS, the crGW-EOS better predicts the pressure of a system, as shown in Fig. 3. There are two reasons for this improvement. First, the low density limits in the crGW-EOS depend solely on the chain density. Second, the high density limits in the crGW-EOS depend largely on the bead density, which includes a corrective factor to account for the effects of bonding. This factor also determines how much 2-bead interaction regions overlap. The effect on the pressure prediction and the appropriate aii for a given density and bond length l can be seen in Fig. 4, which compares the crGW-EOS, GW-EOS, and rGW-EOS. The crGW-EOS accounts for the bond length by adjusting the value of aii so that the chain retains a constant volume regardless of the bond length. Beyond a bond length of 1.0, the beads in the chain no longer overlap. Increasing the bond length has no effect on the volume of the molecule, which we see in Fig. 4 by the saturation of the pressure and aii at bond lengths longer than 1.0. In contrast, both the GW-EOS and rGW-EOS predictions for aii remain completely insensitive to the bonding details of the chains.

FIG. 3.

A comparison of the accuracy of the crGW-EOS for chain molecules with the accuracy of the GW-EOS and the rGW-EOS across the same range of parameters used in Fig. 2. The prediction errors are shown for the crGW-EOS, the rGW-EOS, and the GW-EOS, with dark green circles, magenta pluses, and cyan squares, respectively.

FIG. 3.

A comparison of the accuracy of the crGW-EOS for chain molecules with the accuracy of the GW-EOS and the rGW-EOS across the same range of parameters used in Fig. 2. The prediction errors are shown for the crGW-EOS, the rGW-EOS, and the GW-EOS, with dark green circles, magenta pluses, and cyan squares, respectively.

Close modal
FIG. 4.

A comparison of how each equation of state’s prediction of the pressure (top) fares compared with the simulated pressure of a system. The comparison system consists of a fluid of 3-bead chains at a constant chain density of ρchain = 1.67 with a target pressure of ptarget = 41.89, as a function of bond length. In each plot, the black, red, and green lines depict the results from the crGW-EOS, the rGW-EOS, and the GW-EOS, respectively.

FIG. 4.

A comparison of how each equation of state’s prediction of the pressure (top) fares compared with the simulated pressure of a system. The comparison system consists of a fluid of 3-bead chains at a constant chain density of ρchain = 1.67 with a target pressure of ptarget = 41.89, as a function of bond length. In each plot, the black, red, and green lines depict the results from the crGW-EOS, the rGW-EOS, and the GW-EOS, respectively.

Close modal

To describe binary mixtures composed of homo-oligomer chains, we must develop equations to parametrize the cross-interaction parameter, aij, between two different chains. We first look at athermal mixtures, which are mixtures that have a Flory-Huggins χ-parameter equal to 0. While they are not useful for describing phase equilibria properties, athermal mixtures are useful to verify the cross-interaction parameters from thermal mixtures because in the limit that χ = 0, we should recover the athermal mixture’s cross-interaction parameters. These mixtures behave similarly to ideal mixtures; namely, their enthalpy of mixing is not in excess and the two molecules will have similar chemistry. But, in contrast to ideal mixtures, which upon mixing require constant volume and enthalpy and whose changes in entropy are given only by a combinatorial prefactor, athermal mixtures require only constant enthalpy upon mixing. By considering the free energy of mixing, Travis et al. found a relationship between the DPD parameters and cohesive energy densities.30 This relationship is expressed in terms of the quadratic prefactor and pure-phase density in the DPD equation of state and is given by30 

δiδj2=rcut4αii(ρi0)2aii+αjj(ρj0)2ajj2αij(ρi0)(ρj0)aij.
(25)

This equality describes athermal mixtures when the cohesive energy densities, δi and δj, are equal.

Travis et al. found that when the quadratic cross prefactor, αij, is chosen to be equal to αiiαjj, the cross-interaction parameter can be computed as30 

aij=(ρi0)2aii+(ρj0)2ajj2(ρi0)(ρj0).
(26)

In the development of the rGW-EOS, the quadratic prefactor, α, was replaced by a function, h(a, ρ). This function contains all the terms that are multiplied by 2 in the rGW-EOS.31 We modify the crGW-EOS in the same way and extract the quadratic prefactor for homo-oligomer chains, which is given by

h(a,l,m,ρchain0)=FB2ρchain0B2chainaii,m,la+Fαmρchain0(mρchain0)2×1ηam1m+βlnyRlmρchain0.
(27)

By substituting h(a,l,m,ρchain0) into the rGW-EOS cross-interaction parameter equation, we find the cross-interaction parameter between beads in unlike homo-oligomer chains in the athermal limit by reference to each chain’s pure-phase density properties.

To compare the effectiveness and validity of our new cross-interaction equation against Eq. (26), the GW-EOS, and the rGW-EOS, we ran NpT and NVT simulations of the following fictitious mixtures: Vh1Vh2, Vh1Vo2, and Vh2Vo2. These mixtures and the simulation methods are described in Sec. II C. If the equations of state are robust, the Vh1Vh2 mixture should exhibit ideal thermodynamic mixing behaviors in the athermal mixing limit since both chains represent the same molecule type. Figure 5 demonstrates these behaviors for the Vh1Vh2 mixture, and Fig. 6 demonstrates these behaviors for the Vh1Vo2 and Vh2Vo2 mixtures (i.e., the excess enthalpy of mixing, excess volume of mixing, and excess internal energy of mixing show very small deviations from ideal behavior). The crGW-EOS shows the smallest deviation from ideal behavior, but rGW-EOS and GW-EOS also come close to reproducing athermal mixing behavior. However, where rGW-EOS and GW-EOS deviate significantly is in NVT simulation results, shown in the top left plot of Fig. 5, where pressures show up to an 11% relative error from our target. In contrast, the crGW-EOS-parametrized NVT simulation results for pressure show up to a 2.5% relative error from our target. Comparing Figs. 6 and 5 illuminates an argument for the mixture of Vh1Vo2 that is nearly identical to that of Vh1Vh2. Small deviations from ideal mixture are again observed for each quantity considered. Furthermore, the crGW-EOS yields the smallest deviations from ideal behavior in enthalpy, while rGW-EOS and GW-EOS still diverge far from the correct state point. Figure 6 shows the excess thermodynamic properties of the Vh2Vo2 mixture compared to the Vh1Vo2 mixture. This comparison shows that both mixtures produce small, similar deviations from the expected state point of the system. Therefore, our results indicate that the crGW-EOS is producing nearly athermal behavior based on its small deviations from the desired state point and ideal mixing in nearly all cases. Ultimately, Figs. 5 and 6 ensure that the crGW-EOS and Eq. (27) can parametrize different molecular representations properly.

FIG. 5.

Thermophysical properties for the Vh1Vh2 binary mixture as a function of the mole fraction of Vh1. The properties studied are the relative error in pressure (top left), relative excess enthalpy (bottom left), relative excess internal energy (top right), and relative excess molar volume (bottom right). Solid lines and dotted lines depict results from NpT and NVT simulations, respectively. The magenta, cyan, and purple correspond to the crGW-EOS, rGW-EOS, and GW-EOS, respectively.

FIG. 5.

Thermophysical properties for the Vh1Vh2 binary mixture as a function of the mole fraction of Vh1. The properties studied are the relative error in pressure (top left), relative excess enthalpy (bottom left), relative excess internal energy (top right), and relative excess molar volume (bottom right). Solid lines and dotted lines depict results from NpT and NVT simulations, respectively. The magenta, cyan, and purple correspond to the crGW-EOS, rGW-EOS, and GW-EOS, respectively.

Close modal
FIG. 6.

Thermophysical properties for two binary mixtures, Vh1Vo2 and Vh2Vo2, as a function of the mole fraction of either Vh1 for the Vh1Vo2 mixture (left) or Vh2 for the Vh2Vo2 mixture (right). The properties studied are the relative error in pressure (top), relative excess enthalpy (middle top), relative excess internal energy (middle bottom), and relative excess molar volume (bottom). Solid lines and dotted lines depict results from NpT and NVT simulations, respectively. The magenta, cyan, and purple correspond to the crGW-EOS, rGW-EOS, and GW-EOS, respectively.

FIG. 6.

Thermophysical properties for two binary mixtures, Vh1Vo2 and Vh2Vo2, as a function of the mole fraction of either Vh1 for the Vh1Vo2 mixture (left) or Vh2 for the Vh2Vo2 mixture (right). The properties studied are the relative error in pressure (top), relative excess enthalpy (middle top), relative excess internal energy (middle bottom), and relative excess molar volume (bottom). Solid lines and dotted lines depict results from NpT and NVT simulations, respectively. The magenta, cyan, and purple correspond to the crGW-EOS, rGW-EOS, and GW-EOS, respectively.

Close modal

Finally, Figs. 5 and 6 show that for all three equations of state, the excess enthalpy of mixing, excess volume of mixing, and excess internal energy of mixing deviate from ideal behavior. But, while none of the parametrizations produce perfectly ideal mixing behavior, they produce appropriate athermal mixing behavior. Systems parametrized with GW-EOS and rGW-EOS are consistently incorrect at a nearly constant state point, even though they are at the wrong state point. This indicates that the large excess thermodynamics of mixing for rGW-EOS and GW-EOS observed in the Vh1Vo2 mixture were a result of the large change in the state point. Around that same state point, the crGW-EOS shows similar fluctuations. This demonstrates that a consistent state point is crucial to properly recover athermal mixture parameters. Overall, in properly parametrizing an athermal mixture, the crGW-EOS outperforms the GW-EOS and rGW-EOS methods because it not only produces similar deviations from ideal behavior but also maintains a consistent and accurate pressure/density.

An athermal description of a binary mixture cannot accurately describe the true binary mixture partitioning behavior. There are further chemical and size effects to account for when describing the mixing behavior between two unlike molecules. To properly describe the cross-interaction strength between the two different beads, we accounted for the free energy of mixing. This was accomplished by relating the Flory-Huggins χ-parameter to the binary cross-interaction strength, aij. To compute the χ-parameter of a mixture, we extended the approach developed with the rGW-EOS for single bead molecules to properly account for the effects of bonding on the free energies used.31 For mixtures where the particles have different sizes, the χ-parameter may be expressed in terms of transfer free energies as31,46

χvref2RTρrefΔGi,jTvi+ΔGj,iTvj,
(28)

where ΔGi,jT is the free energy of transfer of a pure-phase particle, i, into a fluid comprising pure-phase particles, j, at a temperature, T. vi is the pure-phase volume of particle i. The χ-parameter may then be expressed in terms of the infinite dilution excess chemical potentials and pure-phase densities as31,46

χ=12RTρrefρi0μi,jexaij,ρj0μi,iexaii,ρi0+ρj0μj,iexaji,ρi0μj,jexajj,ρj0,
(29)

where μi,jexaij,ρj0 is the infinite dilution excess chemical potential of a solute particle, i, that is inserted into a fluid of solvent particles, j, at a density, ρj0. The interaction strength between particles i and j is controlled by the cross-interaction parameter, aij. The χ-parameter is dependent on the infinite dilution excess chemical potential, which is itself dependent on the cross-interaction parameter.

To obtain χ, we used an iterative sampling process, in which we computed the terms for the infinite dilution excess chemical potential present in Eq. (29) by adjusting the value of aij until the desired value of χ is reached. For efficiency, it is best to use a functional relation to pre-compute the infinite dilution excess chemical potential in terms of system parameters. This same approach was taken in the development of the rGW-EOS for single bead molecules, where an expression for the infinite dilution free energy of transfer of a solute particle, i, into a solvent fluid of particles, j, with a pure density, ρj0, was given by the following equation:31 

μi,jexaij,ρj0=1expaij0.00323mjρj0+0.00439×50.9mjρj010.790mjρj0.
(30)

We added the number of beads in the solvent chain, mj, to distinguish the bead density and the pure-phase density (this equation is now in terms of the bead density). With a MUPE of 0.3% for a system with a reference pressure of 41.89 across a wide range of densities, this relation works well for mixtures comprising monomeric representations of molecules.

To extend this equation to predict the free energy of transfer for solute molecules made of more than one bead (mi > 1), we first restrict the extended equation to recover the monomeric result in the limit of one bead per solute chain

limmi1μi,jex,(mi,mj)=μi,jex,(1,mj).
(31)

This requirement ensures that the crGW-EOS is consistent with the rGW-EOS. The notation for excess chemical potential has been modified to include the number of beads mi,mj in each chain. As the bond length between the beads of a solute chain increases into a regime where neighboring beads do not overlap, the chain’s free energy decouples from its bonding parameter. The system behaves as if mi monomer beads were being inserted into the monomer fluid instead of into one bonded chain. This is described as

limli>1μi,jex,(mi,mj)aij,li,mi,aii,ρj0,mj=miμi,jex,(1,mj)aij,li>1,ρj0,mj.
(32)

In the opposite limit, for the bond length li → 0, all beads in the chain occupy the same space. From an enthalpic point of view, they behave as a single bead with an interaction strength mi times the monomer interaction strength

limli=0μi,jex,(mi,mj)aij,li,mi,aii,ρj0,mj=μi,jex,(1)miaij,ρj0,mj.
(33)

Unlike in the monomeric case, the chain self-interactions impede the chain from reaching certain geometries. This impediment increases the free energy of transfer from a reference ideal gas phase into a liquid phase. In a typical DPD simulation, the effects of chain self-interaction on the free energy of transfer are small compared to those of bonding. The limiting behaviors of the free energy of transfer in Eqs. (32) and (33) are independent of the chain self-interactions.

This assumption, along with the stated limits, allows us to introduce a set of switching functions, (fA, fB), which are independent of the chain self-interaction parameter. This allows us to switch smoothly between the two limits given by Eqs. (32) and (33) to calculate the infinite dilution excess chemical potential of a chain fluid transferring into a solvent fluid

μi,jex,(mi,mj)aij,li,aii,ρj0=miμi,jex,(1,mj)aij,ρj0f  A+μi,jex,(1,mj)miaij,ρj0f  B.
(34)

To determine the correct form for these switching functions, we performed a multistep procedure—(1) reduce the parameter space needed to generate and train the fitting functions by using assumptions, (2) run two sets of simulations in the pseudo-osmotic Gibbs ensemble over the remaining parameter space to determine the infinite dilution free energy of transfer of different chain molecules, i (hexamer, tetramer, trimer, and dimer), into a solvent fluid, j, at several solvent densities, ρj0 (outlined in Sec. II D), and (3) use the data from step (2) to infer the form of and train the switching functions to minimize the error in predicted infinite dilution free energy of transfer. We chose which variables were important for determining the form and training the switching functions (fA, fB) in Eq. (34) because the equation varies over a large parameter space. Several variables were eliminated as the physical properties were evaluated. The pure-phase density does not affect the free energy of transfer in the infinite dilution limit, so it was eliminated. Similarly, the solute chain self-interaction constant only indirectly contributes to the free energy of transfer with no effect on the limiting values of the free energy of transfer; therefore, it was not included. For the solute molecule, all other parameters play an important role in the free energy of transfer.

We computed the ratio of excess chemical potentials of a chain and one of the monomer units in that chain. As expected, the ratio approaches the number of beads in the chain as the bond length increases for each homo-oligomer chain (see Fig. S2 of the supplementary material). Additionally, the ratio approaches the value predicted in Eq. (33) in the short bond-length limit. Note that only the solute chain varied in the number of monomers; the solvent molecules were only 1-bead chains. Ultimately, our data allowed us to obtain the switching functions (needed to switch between the stated limits for the bond length in this case), which follow a similar behavior to the switching functions [Eqs. (3) and (4)] used in the rGW-EOS for pressure. Our new switching functions are

f  Ali=c1li21+c2li2,
(35)
f  Bli=11+c3li3.
(36)

To train the equations to the generated data and minimize the maximum relative error, several error minimization methods were used, including a grid scan over initial coefficients. To improve fitting characteristics, only bond lengths ranging from li = 0.25 to li = 0.6 were used. For this training data, we used the fitting coefficients c1 = 1.4851, c2 = 0.46562, and c3 = 2.8239, which yielded a MUPE of 2.094% and a maximum relative error of 3.675%.

As a first test of our new equation for the free energy of transfer, we used the athermal limit of Eq. (29) along with Eq. (34) with the fitted switching functions [Eqs. (35) and (36)] to compute the athermal cross-interaction parameter. We compared the results with the athermal mixing parameter generated by Eqs. (25) and (27). The percent deviation between the two results was computed and then histogrammed (see Fig. S3 of the supplementary material). For this test, we computed every athermal interaction parameter for each combination of pure-phase densities ranging from ρ0 = 1.0 to 9.0 in steps of 1.0, bond lengths ranging from l = 0.2 to 0.7 in steps of 0.05, and number of beads ranging from m = 1.0 to 8.0 with a reference pressure of 41.8, a reference density of ρref. = 5.0, and a reference volume of vref. = 90 Å3, corresponding to the water reference fluid commonly used in DPD simulations. Combinations where the self-interaction constant was not between 5.0 < aii < 60.0 were left out because they likely correspond to non-physical systems. The aij for athermal mixing generated from the thermal mixing equations is, on average, slightly lower than the values generated from the athermal mixing equation. A majority of the aijtherm(χ=0) values were within 2% of the athermal cross-interaction parameter. This suggests that Eq. (34), with the fitted switching functions, provides satisfactory agreement with simulation data over the studied range of parameters at a standard pressure of 41.89 and a density of 5. Note that Eqs. (30), (35), and (36) must be reconstructed for other reference pressures and densities.

From the ratios in excess chemical potential, we tested how the percentage error for each test prediction varies as a function of three different variables: chain length, bond length, and density (see Fig. 7). For this analysis, we chose a single chain length and bond length, then varied the density; we also chose a single bond length and density, then varied the chain length; finally, we also chose a single density and chain length, then varied the bond length.

FIG. 7.

The error in the ratio of the predicted infinite dilution excess chemical potential of the chain to that of a monomer in the chain compared with simulation results as a function of density (red circles), number of beads in the chain (black squares), and bond length (blue diamonds).

FIG. 7.

The error in the ratio of the predicted infinite dilution excess chemical potential of the chain to that of a monomer in the chain compared with simulation results as a function of density (red circles), number of beads in the chain (black squares), and bond length (blue diamonds).

Close modal

Overall, the error found with any of the parameters is within ±4%, so the crGW-EOS is doing well at computing the excess chemical potential. However, the error seems to be trending toward a larger absolute error with longer chains and interestingly, this error goes from a positive value at low values of chain length to a negative error at chain lengths of 4 or more beads. The error does not vary too much with a changing bond length, and it decreases with an increase in density.

To compare the effects that the infinite dilution excess chemical potential, switching functions, and the χ-parameter have on the prediction of the upper-critical solution temperature (UCST) and mixing behavior of binary mixtures, we evaluated the liquid-liquid equilibria (LLE) behavior of three fictional binary mixtures: a Vh1Vh2 mixture, a Vh1Vo2 mixture, and a Vh2Vo2 mixture. For Vh1Vh2 systems, we expected an equimolar critical composition for DPD fluids, and for Vh1Vo2 and Vh2Vo2 systems, we expected a non-equimolar critical composition. DPD conflates particle size and the self-interaction parameters, and this conflation drives the behavior of the mixtures. We expected that if the equation of state adapted well to a different number of beads, then the latter two systems would have nearly identical phase behavior. The parameters used for these simulations are outlined in Sec. II E.

Simulation results are shown in phase diagrams for the Vh1Vo2 and Vh2Vo2 mixtures (see Fig. 8). The results for the Vh1Vh2 mixture can be seen in Fig. S4 of the supplementary material. χ−1 does not differ much between the three equations of state.

FIG. 8.

Composition phase diagram for the Vh1Vo2 mixture (top data) and Vh2Vo2 (bottom data), where the cross-interaction strength was slowly increased until the fluids become totally miscible at χc1 for each equation of state method: crGW-EOS (black), rGW-EOS (red), and GW-EOS (green). The data for the Vh1Vo2 mixture are shifted up by 0.2χ−1 for clarity.

FIG. 8.

Composition phase diagram for the Vh1Vo2 mixture (top data) and Vh2Vo2 (bottom data), where the cross-interaction strength was slowly increased until the fluids become totally miscible at χc1 for each equation of state method: crGW-EOS (black), rGW-EOS (red), and GW-EOS (green). The data for the Vh1Vo2 mixture are shifted up by 0.2χ−1 for clarity.

Close modal

The differentiating behavior between the three equation of state approaches becomes more pronounced for the systems where each molecule has a different pure-phase density. In these systems, the crGW-EOS has a much higher UCST than the rGW-EOS and the GW-EOS approaches, which largely agree with each other for the Vh2Vo2 and Vh1Vo2 mixtures. While the phase behavior for both the Vh1Vo2 and the Vh2Vo2 mixtures is similar when using the crGW-EOS method, the rGW-EOS and GW-EOS methods approach a much higher UCST for the Vh2Vo2 mixtures than for the Vh1Vo2 mixtures. This demonstrates that the rGW-EOS and the GW-EOS representations fail to remain consistent as the number of beads in the molecule changes, while the crGW-EOS remains consistent even with a change in the number of beads used to represent one of the species. For chain molecules, the crGW-EOS yields better results than the rGW-EOS or the GW-EOS methods.

In this section, we study a real binary mixture made up of homo-oligomer chains to validate our devised methods for real systems. To remain consistent with the rGW-EOS development process, we chose to simulate the liquid-liquid equilibria behavior of mixtures of various simple n-alkanes with acetic anhydride, for which there is a wealth of experimental information in the literature.47 Specifically, to test the strength of the crGW-EOS approach, we will study the mixtures of n-hexane/acetic anhydride and n-octane/acetic anhydride. To study these mixtures, we model each alkane molecule with a dimeric representation (corresponding to either 3 or 4 carbon atoms per bead) and the acetic anhydride with both monomeric and dimeric representations.

The COSMOTherm software40,41 was used to compute the transfer free energies and volumes for each of the bead types present in the simulations, as described by Liyana-Arachchi et al. and Tang et al.11,31 Bead volumes were estimated for each bead type at a temperature of T = 298 K with the volumes assumed to be independent of temperature. We estimated the free energy of transfer of each bead type at different temperatures by weighing specific portions of the molecule in the COSMOTherm calculations such that only the molecular surface corresponding to the relevant bead in the molecule was considered. Equation (28) is then used to compute the χ-parameter for each bead type which is then used along with the parameterized form of Eq. (34) to estimate the cross-interaction parameter between each bead type considered in each simulation. Self-interaction parameters for each molecule are computed using the pure-phase molecular density for each molecule type along with Eq. (14). Because the value of aii depends on the average bond length and bond strength between beads in the system, and because those bonding details change as a function of temperature, aii is no longer independent of the temperature, unlike in the rGW-EOS and the GW-EOS. Further parameterization and simulation details can be found in Sec. II F. The results of these simulations are shown in Fig. 9.

FIG. 9.

Coexistence curves for the LLE of mixtures of (top) acetic anhydride/n-hexane and (bottom) acetic anhydride/n-octane. The green circles and triangles represent mixtures where the acetic anhydride is in a 1-bead representation and the blue squares and down-triangles represent mixtures with a 2-bead representation of acetic anhydride. The filled black diamonds are COSMOTherm predictions.31 The black line in each plot represents the experimental data for each mixture.47 

FIG. 9.

Coexistence curves for the LLE of mixtures of (top) acetic anhydride/n-hexane and (bottom) acetic anhydride/n-octane. The green circles and triangles represent mixtures where the acetic anhydride is in a 1-bead representation and the blue squares and down-triangles represent mixtures with a 2-bead representation of acetic anhydride. The filled black diamonds are COSMOTherm predictions.31 The black line in each plot represents the experimental data for each mixture.47 

Close modal

Figure 9 shows that for the simulated hexane/acetic anhydride(1-bead) mixtures, the alkane rich liquid phase agrees with the COSMOTherm predictions. But, for the acetic anhydride rich phase, the model tends to under-predict the amount of hexane in the system when compared with the COSMOTherm predictions. This will likely result in this system exhibiting a higher UCST than would be predicted by COSMOTherm. The experimental data shown in this plot follow a similar trend where the acetic anhydride rich phase has far less hexane than would be predicted by COSMOTherm. In contrast, the hexane/acetic anhydride(2-bead) simulations tend to yield a lower UCST. In the alkane rich phase, the model tends to over-predict the amount of alkane when compared with the COSMOTherm predictions. In the acetic anhydride rich phase, more hexane is shown than is predicted with the hexane/acetic anhydride(1-bead) mixture.

For the octane/acetic anhydride(1-bead) mixture, the DPD simulations tend to have less acetic anhydride in the alkane rich phase than is predicted by COSMOTherm. Similar to the hexane/acetic anhydride mixtures in the alkane rich phase, the 2-bead acetic anhydride model permits more acetic anhydride than is predicted by the COSMOTherm and more than the acetic anhydride(1-bead) mixtures at the same temperature. In the acetic anhydride rich phase, the same type of behavior occurs as we observed in the hexane/acetic anhydride mixtures, namely, that the 2-bead acetic anhydride models permit more octane than the 1-bead acetic anhydride models.

While none of the DPD models seem to get the same composition as the COSMOTherm predictions, they all adequately simulate the phase behavior at each temperature with the acetic anhydride(1-bead) models performing the best. Additionally, each model outperforms the simulations done by Liyana-Arachchi et al.,31 despite the crGW-EOS having a larger range of errors in predicted parameters. This outperformance is likely because the multi-bead representations are usually better representations of these long alkanes where shape and steric effects begin to contribute to the phase behavior and interaction energies of these systems.

The models with the 2-bead representations of acetic anhydride do not show the same LLE phase behavior as the 1-bead acetic anhydride mixtures. The likely cause is the crGW-EOS model itself. Inherent in the prediction of each parameter is a certain amount of error. Recall that for self-interaction parameters we already observed that there could be up to an 11.5% deviation in predicted pressure and simulated pressure for any set of data. Furthermore, the equations for the cross-interaction parameters tend to under-predict the strength of the cross-interaction parameter and are likely to have an error of up to 2.5% according to our previous analysis. The fluctuation in LLE behavior between these 2 models could be explained by these errors.

In this work, we extended the developments for the rGW-EOS to account for the pressure reduction when beads are linked together to form a homo-oligomer chain. Using TPT1, we derived and modified an equation. Aided by insights from the creation of the rGW-EOS, we generated the crGW-EOS, given in Eq. (14). The crGW-EOS is able to predict the pressure of a homo-oligomer chain fluid with a MUPE of 2.5% and a relative error absolute values lower than 11.5% over the ranges of system parameters commonly considered in DPD simulations. For parametrizing chain fluids, this marks an improvement of several orders of magnitude over using the GW-EOS or the rGW-EOS. After developing the crGW-EOS, we developed a set of combining rules for describing interaction between two unlike chains by combining Eqs. (25) and (27) for athermal mixtures, and combining Eqs. (29), (34)–(36) for thermal mixtures. Our equations for thermal mixtures have a small tendency to under-predict the cross-interaction parameters, which leads to an under-prediction of the UCST. Nonetheless, for chain fluids, these new developments for the athermal and thermal mixing rules represent an improvement over the methods used for the rGW-EOS and GW-EOS. Finally, to remain consistent with the rGW-EOS, the new equations were developed so that the rGW-EOS version of each equation is recovered in the limit of 1 bead per chain. Unlike the rGW-EOS work, the simulation protocol is not reported here because the protocols developed for the rGW-EOS are sufficient, and the conclusion remains unchanged.31 

Future work should extend this framework to study hetero-oligomer chains and branched chains so a much larger variety of molecules may be easily and accurately parametrized without the need for many expensive top-down calculations.

See supplementary material for additional data, parameters, and figures.

This research was supported by the National Science Foundation (Grant No. CBET-1159837), the Procter & Gamble Company, and the University of Minnesota Disability Resource Center through access assistants for Dr. Minkara; specifically, Tanner Lambson, Natalie Guse, and Connor Venteicher. Computational resources were provided by the Minnesota Supercomputing Institute. Preliminary research on DPD simulations of chain fluids by Thilanga Liyana-Arachchi is also acknowledged. Finally, we thank Michael Schmidt (P&G) for running the COSMOTherm calculations used to parameterize the n-alkane/acetic anhydride simulations.

1.
P. J.
Hoogerbrugge
and
J. M. V. A.
Koelman
,
Europhys. Lett.
19
,
155
(
1992
).
2.
J. M. V. A.
Koelman
and
P. J.
Hoogerbrugge
,
Europhys. Lett.
21
,
363
(
1993
).
3.
S. J.
Marrink
,
H. J.
Risselada
,
S.
Yefimov
,
D. P.
Tieleman
, and
A. H.
de Vries
,
J. Phys. Chem. B
111
,
7812
(
2007
).
4.
W.
Shinoda
,
R.
DeVane
, and
M. L.
Klein
,
Mol. Simul.
33
,
27
(
2007
).
5.
E. A.
Müller
and
G.
Jackson
,
Annu. Rev. Chem. Biomol. Eng.
5
,
405
(
2014
).
6.
K. A.
Maerzke
and
J. I.
Siepmann
,
J. Phys. Chem. B
115
,
3452
(
2011
).
7.
R. D.
Groot
and
P. B.
Warren
,
J. Chem. Phys.
107
,
4423
(
1997
).
8.
P.
Español
and
P.
Warren
,
Europhys. Lett.
30
,
191
(
1995
).
9.
G.
Kacar
,
E. A. J. F.
Peters
, and
G.
de With
,
Europhys. Lett.
102
,
40009
(
2013
).
10.
A.
Maiti
and
S.
McGrother
,
J. Chem. Phys.
120
,
1594
(
2004
).
11.
X.
Tang
,
W.
Zou
,
P. H.
Koenig
,
S. D.
McConaughy
,
M. R.
Weaver
,
D. M.
Eike
,
M. J.
Schmidt
, and
R. G.
Larson
,
J. Phys. Chem. B
121
,
2468
(
2017
).
12.
A.
Vishnyakov
,
M.
Lee
, and
A. V.
Neimark
,
J. Phys. Chem. Lett.
4
,
797
(
2013
).
13.
M.
Lee
,
A.
Vishnyakov
, and
A. V.
Neimark
,
J. Phys. Chem. B
117
,
10304
(
2013
).
14.
J. C.
Shillcock
and
R.
Lipowsky
,
J. Chem. Phys.
117
,
5048
(
2002
).
15.
M.
Lee
,
R.
Mao
,
A.
Vishnyakov
, and
A. V.
Neimark
,
J. Phys. Chem. B
120
,
4980
(
2016
).
16.
M. A.
Johnston
,
W. C.
Swope
,
K. E.
Jordan
,
P. B.
Warren
,
M. G.
Noro
,
D. J.
Bray
, and
R. L.
Anderson
,
J. Phys. Chem. B
120
,
6337
(
2016
).
17.
M.
Lee
,
A.
Vishnyakov
, and
A. V.
Neimark
,
J. Chem. Phys.
144
,
014902
(
2016
).
18.
R.
Mao
,
M.
Lee
,
A.
Vishnyakov
, and
A. V.
Neimark
,
J. Phys. Chem. B
119
,
11673
(
2015
).
19.
R. D.
Groot
and
T. J.
Madden
,
J. Chem. Phys.
108
,
8713
(
1998
).
20.
R. D.
Groot
,
Langmuir
16
,
7493
(
2000
).
21.
R. D.
Groot
and
K. L.
Rabone
,
Biophys. J.
81
,
725
(
2001
).
22.
K.
Shi
,
C.
Lian
,
Z.
Bai
,
S.
Zhao
, and
H.
Liu
,
Chem. Eng. Sci.
122
,
185
(
2015
).
23.
L.
Rekvig
,
B.
Hafskjøld
, and
B.
Smit
,
Phys. Rev. Lett.
92
,
116101
(
2004
).
24.
K. E.
Novik
and
P. V.
Coveney
,
Phys. Rev. E
61
,
435
(
2000
).
25.
R. D.
Groot
,
T. J.
Madden
, and
D. J.
Tildesley
,
J. Chem. Phys.
110
,
9739
(
1999
).
27.
J. T.
Padding
and
W. J.
Briels
,
J. Chem. Phys.
117
,
925
(
2002
).
28.
N. A.
Spenley
,
Europhys. Lett.
49
,
534
(
2000
).
29.
D. A.
Fedosov
,
W.
Pan
,
B.
Caswell
,
G.
Gompper
, and
G. E.
Karniadakis
,
Proc. Natl. Acad. Sci. U. S. A.
108
,
11772
(
2011
).
30.
K. P.
Travis
,
M.
Bankhead
,
K.
Good
, and
S. L.
Owens
,
J. Chem. Phys.
127
,
014109
(
2007
).
31.
T. P.
Liyana-Arachchi
,
S. N.
Jamadagni
,
D.
Eike
,
P. H.
Koenig
, and
J. I.
Siepmann
,
J. Chem. Phys.
142
,
044902
(
2015
).
32.
M. S.
Wertheim
,
J. Chem. Phys.
85
,
2929
(
1986
).
33.
M. S.
Wertheim
,
J. Chem. Phys.
87
,
7323
(
1987
).
34.
W. G.
Chapman
,
J. Chem. Phys.
93
,
4299
(
1990
).
35.
J. K.
Johnson
,
E. A.
Mueller
, and
K. E.
Gubbins
,
J. Phys. Chem.
98
,
6413
(
1994
).
36.
Q. P.
Chen
,
B.
Xue
, and
J. I.
Siepmann
,
J. Chem. Theory Comput.
13
,
1556
(
2017
).
37.
M. G.
Martin
and
J. I.
Siepmann
,
J. Phys. Chem. B
102
,
2569
(
1998
).
38.
D. A.
McQuarrie
,
Statistical Mechanics
(
University Science Books
,
2000
).
39.
V. I.
Harismiadis
and
I.
Szleifer
,
Mol. Phys.
81
,
851
(
1994
).
40.
A.
Klamt
,
F.
Eckert
, and
W.
Arlt
,
Annu. Rev. Chem. Biomol. Eng.
1
,
101
(
2010
).
41.
A.
Klamt
,
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
1
,
699
(
2011
).
42.
M. G.
Martin
and
J. I.
Siepmann
,
Theor. Chem. Acc.
99
,
347
(
1998
).
43.
A. Z.
Panagiotopoulos
,
N.
Quirke
,
M.
Stapleton
, and
D. J.
Tildesley
,
Mol. Phys.
63
,
527
(
1988
).
44.
H.
Sun
,
J. Phys. Chem. B
102
,
7338
(
1998
).
45.
P. B.
Warren
,
A.
Vlasov
,
L.
Anton
, and
A. J.
Masters
,
J. Chem. Phys.
138
,
204907
(
2013
).
46.
C. M.
Wijmans
,
B.
Smit
, and
R. D.
Groot
,
J. Chem. Phys.
114
,
7644
(
2001
).
47.
I.
Garca de la Fuente
,
M.
Aboy
,
S.
Villa
,
N.
Riesco
,
J. A.
Gonzlez
, and
J. C.
Cobos
,
J. Chem. Eng. Data
47
,
950
(
2002
).

Supplementary Material