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 *B*_{2} 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.

## I. INTRODUCTION

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-grain^{6} 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, *a*_{ij}, to the thermodynamic compressibility in systems comprising particles that interact with the DPD potential.^{7} The GW-EOS is given by

where *p*_{GW} is the pressure of the system, *T* is the absolute temperature, *k*_{B} 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 *a*_{ii} 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 *a*_{ii} 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 *B*_{2} virial coefficient and a set of switching functions to restore virial characteristics expected at low densities^{31} and to adjust how much the quadratic terms contribute to the pressure in the equation of state. The rGW-EOS is given by^{31}

with

and

where *B*_{2}(*a*_{ii}) is the second virial coefficient of a unary DPD fluid with the repulsive force constant of $aii$, *c*_{1} and *c*_{2} 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}

## II. SIMULATION METHODS

### A. DPD potential

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

where *a*_{ij} describes the strength of the interaction between two particles, *i* and *j*, which are separated by a distance, *r*_{ij}. This potential is truncated at a distance in dimensionless DPD units, *r*_{cut}, 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, *k*_{spring}, and a corresponding bond length, *l* = *r*_{bond}/*r*_{cut}. The strength of the bonds does not affect the system pressure if the bonds have a rigidity *k*_{spring} > 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 *k*d-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.

### B. Simulations for developing and validating the crGW-EOS

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, *a*_{ii} (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}

where *U*_{inter} is the potential between two molecules whose centers of mass are a distance, *r*_{12}, apart. The $\u27e8\u2026\u27e9\alpha 1,\alpha 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 *r*_{12} values in step sizes of 0.05 Å/*r*_{cut}. We run these two-box simulations using the same range of system parameters (*m*, *l*, and *a*_{ii}) as with the one-box tuning simulations. The results are averaged across eight independent simulations.

### C. Simulations for athermal mixtures

To validate a cross-interaction term, *a*_{ij}, 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, *V*_{h1}, with a molecular volume similar to that of *n*-hexane (220 Å^{3}) in a mixture with a 2-bead fluid, *V*_{h2}, of the same volume (220 Å^{3}); *V*_{h1} in a mixture with a 2-bead fluid, *V*_{o2}, with a molecular volume similar to that of *n*-octane (275 Å^{3}); and *V*_{h2} in a mixture with *V*_{o2} (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 *V*_{h1} and the volume of a single bead in *V*_{h2}, 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, *r*_{cut}, of 9.086 and a reference pressure, *p*_{ref}, 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 *V*_{h2} and *V*_{o2} 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 *x*_{i} = 0.0 to *x*_{i} = 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.

. | a_{ii}
. | . | |
---|---|---|---|

Model . | V_{h1}
. | V_{h2}
. | a_{ij}
. |

crGW-EOS | 35.0 | 8.57 | 17.3 |

rGW-EOS | 35.0 | 7.43 | 16.1 |

GW-EOS | 32.3 | 7.36 | 15.4 |

V_{h1} | V_{o2} | ||

crGW-EOS | 35.0 | 13.2 | 21.5 |

rGW-EOS | 35.0 | 11.9 | 20.4 |

GW-EOS | 32.3 | 11.6 | 19.4 |

V_{h2} | V_{o2} | ||

crGW-EOS | 8.57 | 13.2 | 10.6 |

rGW-EOS | 7.46 | 11.9 | 9.4 |

GW-EOS | 7.36 | 11.6 | 9.24 |

. | a_{ii}
. | . | |
---|---|---|---|

Model . | V_{h1}
. | V_{h2}
. | a_{ij}
. |

crGW-EOS | 35.0 | 8.57 | 17.3 |

rGW-EOS | 35.0 | 7.43 | 16.1 |

GW-EOS | 32.3 | 7.36 | 15.4 |

V_{h1} | V_{o2} | ||

crGW-EOS | 35.0 | 13.2 | 21.5 |

rGW-EOS | 35.0 | 11.9 | 20.4 |

GW-EOS | 32.3 | 11.6 | 19.4 |

V_{h2} | V_{o2} | ||

crGW-EOS | 8.57 | 13.2 | 10.6 |

rGW-EOS | 7.46 | 11.9 | 9.4 |

GW-EOS | 7.36 | 11.6 | 9.24 |

### D. Simulations for thermal mixtures

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}

where $\rho liq,ji$ and $\rho 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.

Solvent parameters . | ||
---|---|---|

Input density . | Simulated density . | a_{jj}
. |

2.0 | 1.96_{1} | 126.6 |

2.5 | 2.44_{2} | 73.25 |

3.0 | 2.96_{2} | 47.72 |

3.5 | 3.481_{3} | 33.52 |

4.0 | 3.99_{4} | 24.80 |

4.5 | 4.50_{5} | 19.05 |

5.0 | 5.01_{6} | 15.07 |

5.5 | 5.51_{6} | 12.19 |

Solvent parameters . | ||
---|---|---|

Input density . | Simulated density . | a_{jj}
. |

2.0 | 1.96_{1} | 126.6 |

2.5 | 2.44_{2} | 73.25 |

3.0 | 2.96_{2} | 47.72 |

3.5 | 3.481_{3} | 33.52 |

4.0 | 3.99_{4} | 24.80 |

4.5 | 4.50_{5} | 19.05 |

5.0 | 5.01_{6} | 15.07 |

5.5 | 5.51_{6} | 12.19 |

### E. Simulations for miscibility and UCST

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: *V*_{h1} with *V*_{h2}, *V*_{h1} with *V*_{o2}, and *V*_{h2} with *V*_{o2}. 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 *r*_{cut} 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.

Model . | m-beads
. | a_{ii}
. | l_{i}
. |
---|---|---|---|

crGW (hexane) | 1 | 35.0 | … |

2 | 8.57 | 0.39 | |

crGW (octane) | 2 | 13.2 | 0.5 |

rGW (hexane) | 1 | 35.0 | … |

2 | 7.43 | 0.39 | |

rGW (octane) | 2 | 11.9 | 0.5 |

GW (hexane) | 1 | 32.3 | … |

2 | 7.36 | 0.39 | |

GW (octane) | 2 | 11.6 | 0.5 |

Pure-phase density (hexane) | 3.432 | ||

Pure-phase density (octane) | 2.780 |

Model . | m-beads
. | a_{ii}
. | l_{i}
. |
---|---|---|---|

crGW (hexane) | 1 | 35.0 | … |

2 | 8.57 | 0.39 | |

crGW (octane) | 2 | 13.2 | 0.5 |

rGW (hexane) | 1 | 35.0 | … |

2 | 7.43 | 0.39 | |

rGW (octane) | 2 | 11.9 | 0.5 |

GW (hexane) | 1 | 32.3 | … |

2 | 7.36 | 0.39 | |

GW (octane) | 2 | 11.6 | 0.5 |

Pure-phase density (hexane) | 3.432 | ||

Pure-phase density (octane) | 2.780 |

### F. Simulations for acetic anhydride/*n*-alkane mixtures

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 *r*_{cut} = 7.66 Å, and *a*_{ii,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.

### G. Calculation of cavity correlation function

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

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 solver^{45} 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.

## III. RESULTS AND DISCUSSION

### A. Development of the crGW-EOS

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 *B*_{2} 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:

where *p*_{R} 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., *ρ*_{R}*k*_{B}*T*$vs.$ (*ρ*_{R}/*m*)*k*_{B}*T*) 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

##### a. Physical interpretation of the *y*_{R}(*r*) term.

*y*_{R}(*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, *y*_{R}(*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 $[\u2009y2(l)]m\u22121$ [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 $\rho R2\u2202\u2009\u2202\rho Rln[\u2009yR(l)]$ term in Eq. (11) can be as follows:

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/*V*_{R}, we get $\Delta P=\rho R2\u2202A/\u2202\rho R=\u2212\rho R2\u2202wR(l)/\u2202\rho 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 $\Delta P\u223cm\u22121m\rho R2\u2202wR(l)/\u2202\rho R$.

#### 2. Incorporating the low density limit and the role of *B*_{2}

In the low density limit, we expect the pressure of the chain fluid to behave as $P/kBT=\rho chain+B2\rho chain2$, where *B*_{2} 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

##### 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

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

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

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

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\pi 3$ for a large *a*_{ii}, and it is valid for any positive valued *a*_{ii}. This is important because most of the *a*_{ii}-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 *a*_{ii} less than 50. At *a*_{ii} > 50, the two relations diverge.

To find a functional form for *g*(*a*_{ii}, *l*, *m*) and $f\u2009\u2009m,l$, we had to make them independent of one another. To do this, we found a functional relation for $f\u2009\u2009m,l$ because Eq. (16) implied

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

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*(*a*_{ii}, *l*, *m*) by factoring out the contribution from the monomeric portion of the equation

The resultant quantity could be described by the approximation

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),

### B. Fitting and validating the crGW-EOS

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.

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 *a*_{ii} 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 *a*_{ii} 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 *a*_{ii} at bond lengths longer than 1.0. In contrast, both the GW-EOS and rGW-EOS predictions for *a*_{ii} remain completely insensitive to the bonding details of the chains.

### C. Athermal mixtures

To describe binary mixtures composed of homo-oligomer chains, we must develop equations to parametrize the cross-interaction parameter, *a*_{ij}, 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 by^{30}

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 $\alpha ii\alpha jj$, the cross-interaction parameter can be computed as^{30}

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 *aρ*^{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

By substituting $h(a,l,m,\rho 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: *V*_{h1}–*V*_{h2}, *V*_{h1}–*V*_{o2}, and *V*_{h2}–*V*_{o2}. These mixtures and the simulation methods are described in Sec. II C. If the equations of state are robust, the *V*_{h1}–*V*_{h2} 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 *V*_{h1}–*V*_{h2} mixture, and Fig. 6 demonstrates these behaviors for the *V*_{h1}–*V*_{o2} and *V*_{h2}–*V*_{o2} 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 *V*_{h1}–*V*_{o2} that is nearly identical to that of *V*_{h1}–*V*_{h2}. 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 *V*_{h2}–*V*_{o2} mixture compared to the *V*_{h1}–*V*_{o2} 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.

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 *V*_{h1}–*V*_{o2} 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.

### D. Thermal mixtures

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, *a*_{ij}. 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 as^{31,46}

where $\Delta 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 as^{31,46}

where $\mu i,jexaij,\rho 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, $\rho j0$. The interaction strength between particles *i* and *j* is controlled by the cross-interaction parameter, *a*_{ij}. 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 *a*_{ij} 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, $\rho j0$, was given by the following equation:^{31}

We added the number of beads in the solvent chain, *m*_{j}, 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 (*m*_{i} > 1), we first restrict the extended equation to recover the monomeric result in the limit of one bead per solute chain

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 *m*_{i} monomer beads were being inserted into the monomer fluid instead of into one bonded chain. This is described as

In the opposite limit, for the bond length *l*_{i} → 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 *m*_{i} times the monomer interaction strength

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, (*f*_{A}, *f*_{B}), 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

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, $\rho 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 (*f*_{A}, *f*_{B}) 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

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 *l*_{i} = 0.25 to *l*_{i} = 0.6 were used. For this training data, we used the fitting coefficients *c*_{1} = 1.4851, *c*_{2} = 0.46562, and *c*_{3} = 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 < *a*_{ii} < 60.0 were left out because they likely correspond to non-physical systems. The *a*_{ij} 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(\chi =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.

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.

### E. Miscibility and UCST

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 *V*_{h1}–*V*_{h2} mixture, a *V*_{h1}–*V*_{o2} mixture, and a *V*_{h2}–*V*_{o2} mixture. For *V*_{h1}–*V*_{h2} systems, we expected an equimolar critical composition for DPD fluids, and for *V*_{h1}–*V*_{o2} and *V*_{h2}–*V*_{o2} 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 *V*_{h1}–*V*_{o2} and *V*_{h2}–*V*_{o2} mixtures (see Fig. 8). The results for the *V*_{h1}–*V*_{h2} mixture can be seen in Fig. S4 of the supplementary material. χ^{−1} does not differ much between the three equations of state.

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 *V*_{h2}–*V*_{o2} and *V*_{h1}–*V*_{o2} mixtures. While the phase behavior for both the *V*_{h1}–*V*_{o2} and the *V*_{h2}–*V*_{o2} mixtures is similar when using the crGW-EOS method, the rGW-EOS and GW-EOS methods approach a much higher UCST for the *V*_{h2}–*V*_{o2} mixtures than for the *V*_{h1}–*V*_{o2} 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.

### F. LLE of various acetic anhydride/*n*-alkane mixtures

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 software^{40,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 *a*_{ii} 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, *a*_{ii} 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.

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.

## IV. CONCLUSIONS

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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.