Group-contribution coarse-grained molecular simulations of polystyrene melts and polystyrene solutions in alkanes using the SAFT-γ force field

A coarse-grained (CG) model for atactic polystyrene is presented and studied with classical molecular-dynamics simulations. The interactions between the CG segments are described by Mie potentials, with parameters obtained from a top-down approach using the SAFT-γ methodology. The model is developed by taking a CG model for linear-chain-like backbones with parameters corresponding to those of an alkane model, and decorating it with side branches with parameters from a model of toluene, which incorporate an “aromatic-like” nature. The model is validated by comparison with the properties of monodisperse melts, including the effect of temperature and pressure on density, as well as structural properties (the radius of gyration and end-to-end distance as functions of chain length). The model is employed within large-scale simulations that describe the temperature-composition fluid-phase behaviour of binary mixtures of polystyrene in n-hexane and n-heptane. A single temperature-independent unlike interaction energy parameter is employed for each solvent to reproduce experimental


■ INTRODUCTION
As a consequence of its relatively simple chemical structure, polystyrene (PS) has become a ubiquitous model system for fundamental studies aimed at understanding the behavior of more-complex macromolecules and supramolecular assemblies. From the point of view of molecular simulation, PS is a useful benchmark, as it has been studied with fully atomistic, unitedatom, and coarse-grained representations as well as multiscale approaches. From the considerable literature pertaining to fully atomistic or united-atom studies, one can highlight refs 1−7 as exemplifying the varied properties and phenomena studied: structure and dynamics; polymer networks and their gas solubilities and diffusivities; stress response; and the elastic behavior of thin films. In spite of impressive advances in computer power, the computational expense incurred by retaining such a detailed description still limits the system size of the computational cell, which scales exponentially with the molecular weight of the polymer. This becomes relevant as some macroscopic properties (such as phase equilibria) require particularly large systems to properly explore the several coexisting bulk phases. This difficulty can be alleviated through coarse graining; the atomistic detail and the degrees of freedom of the model are reduced by representing groups of atoms as simple segments or "superatoms", where the interactions are modeled with effective force fields. 8 In the current work we focus on modeling PS melts and PS in solution with alkanes, which correspond to scales of length and time where coarsegrained (CG) models provide an appropriate description.
Coarse graining the polystyrene molecule requires a grouping of the force-field centers in such a way that the intermolecular forces (which characterize the overall energetics) and the shape and morphology (which characterize the entropic contributions of the chain molecules) are both duly acknowledged. It is ultimately the balance between energy and entropy, i.e., the free energy, that determines the macroscopic thermophysical quantities. It is in this sense that the SAFT coarse-graining approach excels; the overall shape and connectivity of the original molecules are retained while the free energy is faithfully represented through a direct link to an equation of state.
A comprehensive review of recent CG models of PS has been given by Karimi-Varzaneh et al. 9 They summarize the general characteristics of the various CG models, the methodologies used to obtain the force fields, and comparisons of the performance of each model in terms of the representation of a variety of structural and dynamical properties, of melts (both without deformation and under shear) and solutions. 10−26 An outstanding problematic issue with existing CG models is the lack of representability; the performance of the model deteriorates when it is used for thermodynamic conditions which are different to those pertaining to the states considered in the development of the model parameters (cf. the specific examples discussed in refs 9, 21, 25, and 27).
In contrast to the CG approaches reviewed in ref 9, here we adopt the SAFT-γ top-down methodology reviewed in ref 8. An analytical equation of state (the most recent version 28 of the statistical associating fluid theory, SAFT 29,30 ) is used to develop molecular-force-field parameters from macroscopic thermodynamic information, both for the solvents and the selected molecules used to represent the PS moieties. Under the SAFTγ description, the molecules are modeled as chains of spherical segments characterized with an isotropic intermolecular potential of the Mie form. While any engineering equation of state (EoS) can be used as a tool for correlating and extrapolating macroscopic thermophysical data, SAFT is one of the very few rigorous molecular-based EoSs that is sufficiently accurate to represent the intermolecular potential on which the corresponding molecular model is based.
A number of versions of the SAFT EoS have been developed, distinguished primarily by the interaction potential employed to describe the intermolecular forces between the unbonded segments (see, e.g., ref 31 for a review). In our current work we employ a recently developed group-contribution version, referred to as SAFT-γ Mie or, simply, SAFT-γ, 28 which is based on the parent SAFT-VR Mie 32 formulation and employs the Mie potential to represent the intersegment dispersion interactions. An abridged description of SAFT-VR Mie is presented in refs 33 and 34. A unique feature of the analytical expressions contained within both SAFT-VR Mie and SAFT-γ is the accuracy in linking macroscopic thermodynamic properties to the underlying molecular force fields. It follows that SAFT-γ allows a direct link to be made between the experimental thermophysical property data (e.g., densities, vapor pressures, etc.) and the parameters of a force field that corresponds to an effective pairwise potential exhibiting the same macroscopic behavior when employed in molecular simulation.
Modeling the phase behavior of polymer solutions using molecular simulation is a multifaceted challenge. Models for polymers must be detailed enough to incorporate the correct physics (i.e., connectivity, energetics, etc.) but should not be onerous in terms of computation and should hence incorporate some degree of coarse graining; obtaining the right balance is the crux of the problem. In our current paper we describe the development of a CG model of PS using SAFT-γ and compare the resulting thermophysical and structural properties with the corresponding properties of higher-fidelity models and real PS, as determined experimentally. The fluid-phase behavior of solutions of PS in n-hexane and n-heptane solvents is explored.

■ POLYMER SOLUTIONS
A generic temperature−volume fraction (composition) (T, Φ) phase diagram for polymer solutions 35−37 is provided in Figure  1. The regions highlighted represent liquid−liquid (LL) demixing regions, each bounded by an upper critical solution temperature (UCST) or a lower critical solution temperature (LCST). Typically, if the molecular weight of polymer is increased, the UCST tends to rise and the LCST tends to lower; for higher molecular weights, both regions can merge to form a so-called "hourglass" demixing phase diagram. The appearance of the LCST in the higher temperature region is well understood as an effect due to the compressibility of the solution at such conditions. References 35−38 offer a deeper discussion of polymer−solvent diagrams. Here, we use the SAFT-γ CG models for PS and solvents to explore this type of fluid-phase behavior via molecular dynamics (MD) simulation. We note in passing that it is not appropriate to use the SAFT-γ Mie EoS directly to calculate the phase diagrams for branched chains since the details of structure and connectivity of the polymer are not explicitly taken into account.
The experimental values of the UCST and LCST for the systems of interest 35,39 are summarized in Table 1; the experiments were carried out for volume fractions of polymer between 0 and 0.4 in the case of polymers with a narrow molecular weight distribution (M w /M n < 1.06, where M n is the number-average molecular weight). We adopt the nomenclature "PSx" to refer to polystyrene with a molecular weight M w of x g mol −1 . The fluid-phase behavior of solutions of PS4800 in n-hexane is characterized by an "hourglass" diagram, wherein demixing appears in the interval (0.07 ≤ Φ ≤ 0.32) for all temperatures. 39 The other systems exhibit upper and lower The demixing regions become larger with increasing polymer molecular weight; the UCST increases and the LCST decreases (as indicated by the arrows) until both regions merge into an "hourglass"shaped diagram (represented by the gray shading).

■ MODELS
The molecular models adopted here are coarse-grained chains of tangent spherical segments, interacting via a Mie potential. We use the SAFT-γ Mie EoS 28,32 to obtain the intersegment potential parameters relating to the CG segments. Each molecule is represented by n spherical segments characterized by a size σ, closely related to an average diameter. We model the intersegment interactions with a particular case (λ, 6) of the Mie potential: where r is the distance between centers of a pair of segments, ϵ is the depth of the potential energy well, λ is the repulsive exponent, and is a constant defined (for the (λ, 6) Mie) as The adoption of the particular Mie (λ, 6) form (in preference to the more generic (λ r , λ a ) form) follows from the demonstration 40 that for the description of fluid-phase behavior the repulsive and attractive exponents are not completely independent; this inter-relationship reduces the degrees of freedom, thereby allowing the attractive exponent to be fixed to the London dispersion 41 value of six. The potential parameters, σ, ϵ, and λ, are obtained in one of two ways, depending on the availability of suitable experimental data. In the first instance, the parameters are estimated in order to best represent the experimental vapor−liquid phase-equilibrium data, as described in refs 8 and 42. Alternatively, the parameters can be determined in expedited fashion from corresponding-states correlations, linking them to the experimental acentric factor, critical temperature, and liquid density, respectively; the details of the latter methodology are given in refs 43−45.
Solvents. Two alkane solvents are examined in our current study: n-hexane and n-heptane. The solvents are modeled as homonuclear chains comprising n = 2 segments, with each segment chosen to represent a group of three to four carbon backbone atoms. The intersegment parameters are obtained following the procedure set out in refs 43 and 44 and are given in Table 2. The same set of parameters is used for the SAFT-γ theoretical phase diagrams as for the MD simulations presented in the Results section.
Polystyrene. In group-contribution approaches one can exploit the transferability of CG models for small molecules (or chemical functional groups of atoms) in the development of models of larger or more complicated molecules; 28,46,47 this lies at the heart of the CG philosophy. The SAFT-γ methodology has already been successfully applied in this way to represent the phase equilibria and structure of linear polymers and surfactants. 48−51 Conversely, a large molecule can be subdivided into functional groups chosen to simplify the coarse graining; this is the approach taken here to build our polystyrene CG model. The procedure is summarized schematically in Figure 2. The atomistic polystyrene chain can be seen as an aliphatic backbone decorated with aromatic branches. Our CG "backbone" is a long alkane with parameters estimated the using SAFT-γ Mie methodology using target experimental vapor−liquid equilibria data relating to n-C 9 H 20 , n-C 12 H 26 , and n-C 15 H 20 , as presented in ref 48. The original alkane model comprised two types of segments: "end" segments, representing CH 3 −CH 2 −CH 2 −, and "middle" segments, representing −CH 2 −CH 2 −CH 2 −; intramolecular interactions (bonds and angular restrictions) derived from fitting distribution functions obtained from atomistic simulations were also included. For the sake of simplicity, we reduce the backbone to a homonuclear chain of "middle" segments ( Figure 2b,d) and omit the intramolecular interactions, resulting in a fully flexible chain of tangent segments. The "backbone" alkyl chain parameters are given in Table 2.
The moiety used to represent the aromatic branches is half (one segment) of our (dimer) model of a toluene molecule, as shown in Figure 2c,d; the molecular parameters for this toluene model are obtained following ref 43 and correspond to those of the "branch" aromatic moiety given in Table 2. a Unlike interaction parameters, k ij , segment diameter, σ; energy well depth, ϵ (after applying the corresponding k ij parameter); and the repulsive Mie exponent, λ, respectively. The polymer and solvent molecules are constructed from these building blocks as described in the text.

Article
In summary, the monomer (repeat unit) of our PS model is represented by two segments: one aliphatic in nature and another aromatic in nature, thereby accounting for the important physical features required for a model at this level of coarse-graining. In comparison to fused models, tangential models effectively include more carbons per monomer than is implied by the original molecule or atomic grouping on which the CG model is based (cf. the discussion in ref 47). This is accounted for insofar as the segments of our model represent 3 carbons for the aliphatic segment and 3.5 carbons for the aromatic segment, tallying to 6.5 carbons per monomer; in comparison, there are a total of 8 carbons in the monomer of an actual PS molecule. This apparently odd choice of attaching only half a toluene molecule to every three backbone carbons (represented by a backbone segment) should not be taken literally: top-down CG models do not necessarily reproduce the atomistic detail but incorporate information about the effective average interactions and general shape and connectivity of the molecule. We also note that in our model the aromatic branches are in effect, randomly oriented in relation to the backbone, whereby the model best represents an atactic PS molecule. The parameters relating to the intersegment potentials for the PS model are given in Table 2.
Unlike Interactions for Mixtures. The combining rules used to represent the cross-interactions between unlike segments i and j are (the Lorentz rule) for the diameters, for the repulsive exponent, 32 and (Berthelot-like rule) for the energy, 32 where the adjustable binary-interaction parameter, k ij , is taken as zero unless otherwise stated. Reliable CG models for pure systems do not necessarily lead to a good description of their mixtures, since the standard combining rules (with k ij = 0) correspond to mathematical averages that strictly speaking have a sound physical justification only for simple, small, near-spherical molecules that interact only through London dispersion interactions. 52 For more-complicated molecules, such as those considered in our current study, adjustment of the binary-interaction parameter k ij in eq 5 allows one to account for this extra complexity, providing ultimately for a better description of the experimental behavior. Following the SAFT-γ methodology, the estimation of binary-interaction parameters is typically carried out by minimizing an objective function constructed as the difference between an EoS prediction and experimental information related to fluid-phase equilibria, e.g., isothermal pressure−composition or isobaric temperature−composition data, or to other thermodynamic properties such as volumes of mixing. 37,53−56 However, in the SAFT theory we employ, it is implicit that the chains representing the molecules are linear. In this sense, there is a clear mismatch between the theory and the simulations so that this procedure is not directly appropriate for our current study. Rather, any required adjustment to the k ij values can be judged only from the MD simulations themselves.
In our simulations two unlike interactions between PS and the solvent (n-hexane or n-heptane) are included: backbone− solvent k bbn−sol , and branch−solvent k bch−sol . In practical terms, obtaining both of the corresponding binary-interaction parameters simply by adjusting the MD simulation data to capture the experimental behavior would be prohibitively timeconsuming. However, the influence of these two parameters is partially coupled, and one would expect some degeneracy insofar as different combinations of these two parameters will lead to essentially the same behavior. We therefore assign a value to one unlike-interaction parameter and adjust the other. The backbone should feature an affinity to aliphatic solvents; to accentuate this affinity, we assign a small negative value, k bbn−sol = −0.01, for the backbone−solvent interaction in all our systems. The remaining unlike-interaction parameter, k bch−sol , is adjusted directly in the simulations. In order to bound the k bch−sol value for each solvent, simulations are carried out at two key temperatures: one at which immiscibility is expected, to obtain the lower-bound of k bch−sol , and one where full miscibility is the target; the upper bound is identified as the maximum value of k bch−sol that fulfills this condition. A third simulation at another temperature is used where necessary. The fluid-phase behavior of the mixture is characterized using histograms of the local compositions; these are discussed further in the Results section.
The branch−solvent unlike-interaction parameters used in our current work are k bch−sol = 0.055 when the solvent is n-hexane and k bch−sol = 0.045 when the solvent is n-heptane. The unlike parameters are summarized in Table 2. We emphasize that the unlike-interaction parameters for both solvents are obtained using the PS2030 systems and transferred without further modification to the PS4800 systems for use over a wide range of temperatures and compositions.

■ SIMULATION DETAILS
Three types of systems are simulated: pure solvents, PS melts, and PS solutions. All of the MD simulations are performed using HOOMD-blue 57−59 on single NVIDIA Tesla K40 GPUs. The fluid-phase equilibria and single-phase systems are studied under canonical (NVT) and/or isobaric−isothermal (NPT) conditions. Neighboring segments on the same chain interact with each other only through a harmonic-spring potential characterized by the spring constant k H ; segments separated by more than one bond interact via a Mie potential, as described earlier. In this way our model retains a link with the SAFT description. The spring constant, relating to the bond between adjacent segments, is set to k H = 20 000 kJ mol −1 nm −2 . The time step is Δt = 0.01 ps in all cases, and the cutoff distance is set to 28 Å. The NVT simulation is accomplished using the Martyna−Tobias−Klein (MTK) integration of the equations of motion 60,61 based on the Nose−Hoover thermostat, 62,63 while the NPT simulation is performed using the MTK barostat− thermostat. 60,64,65 Solvent systems with N mol = 5000 molecules (N tot = 10 000 segments) are simulated over t = 80 ns, corresponding to wall-clock times between 7 and 14 h. The most computationally demanding simulations are the mixtures, with around N tot ≈ 10 5 segments. Some of these systems require t = 400 ns of simulation time, corresponding to 800 h of wall-clock time. Specific details for each type of system are given in the Supporting Information.  68,69 Although the explicit inclusion of an interface inevitably requires the use of relatively large system sizes, conventional computational capabilities are sufficient in the case of the system under study here. The use of direct coexistence methods circumvents issues associated with insertion (and deletion) of large particles at liquid-like densities required to achieve equilibrium states through GEMC simulations. In Figure 3, we present the vapor−liquid coexistence envelopes and the vapor-pressure curves obtained for the fluid-phase equilibria of n-hexane and n-heptane by MD simulation. For comparison, we include the smoothed experimental information obtained from NIST 66,67 as well as SAFT-γ 28 EoS calculations, using the same models (Table 2) as in the MD simulations. The average simulated values are obtained using the analysis tools from GROMACS. 70 A good correspondence between theory, simulation, and experiment is apparent from both plots,

Macromolecules
Article although the quality of the descriptions with both the theory and simulation is better in terms of the saturation densities (Figure 3a) than the vapor pressures (Figure 3b), where small deviations of both the theory and simulation from experiment can be seen at higher temperatures.
To assess the compressibility of the CG models in the liquid region, NPT simulations of single-phase systems are performed. Three isotherms are tested for each system starting from pressures close to the corresponding saturation value up to 100 bar. The densities obtained after equilibration are plotted in Figure 4, where the corresponding NIST data 66,67 and SAFT-γ calculations are included for comparison. A good correspondence is observed between all of the data sets, with a slight deterioration for the higher temperature isotherms, which is expected due to the proximity to the critical region.
Polystyrene Melts. To test the structural properties of the PS model, monodisperse systems of PS melts corresponding to different chain lengths are simulated, as detailed in Table 3. The equilibrium density is first obtained from the NPT simulations using a molar mass for the styrene unit of M w = 104.15 g/mol.   the backbone−branch angle, defined as the angle formed by consecutive bbn−bbn−bch triads; and (c) the dihedral angle, defined as the angle between the projections of consecutive bbn−bch pairs on the plane perpendicular to the connecting bbn−bbn diad. Our distributions are obtained from the PS system with N chains = 500 polymers of N mono = 300 styrene units simulated at T = 500 K. For comparison, we include the corresponding results from CG models with similar structure and connectivity: Hsu et al., 21 CG model simulated at a temperature of T = 300 K; Rao et al., 22 CG model for chains with 40 monomers at a temperature of T = 300 K; and Rosch et al., 25 CG ("locally shifted higher resolution") model at a temperature of T = 500 K. All cases are normalized with respect to the corresponding highest peak.

Article
Our simulated densities vary from 900 to 918 kg/m 3 for systems with 10 ≤ N mono ≤ 300, where N mono refers to a monomeric moiety, composed in this case of a dimer of a backbone CG segment and a branch CG segment. The density slightly increases with increasing chain length, in agreement with the trend reported in the literature. 10,17,24 The reported experimental density 15,71 for the 10-mer melt is 895 kg/m 3 , while a density of ∼ 960 kg/m 3 was determined Karimi-Varzaneh et al. 9 for fully atomistic 10-mers, at similar conditions. Once the NPT simulations are equilibrated, NVT simulations are carried out to obtain the ensemble average endto-end distances ⟨R⟩ (calculated from the initial and final backbone segments) and radii of gyration ⟨R g ⟩; these are also reported in Table 3. As a test of the system size, the larger polymer, N mono = 300, is simulated with N chain = 200 and N chain = 500 polymer molecules. The averages obtained for both systems are found to differ by less than 1%.
Plots of the end-to-end distances and radii of gyration as a function of the chain length of the polymer are presented in Figures 5 and 6 together with values for other CG models taken from the literature for comparison. The structural results for our CG polystyrene are quantitatively similar to those obtained   15 Qian et al., 20 Fritz et al. 19 (in agreement with their atomistic simulations, not shown), Hsu et al., 21 and Carbone et al., 27 both for their CG and atomistic models. Note that in the CG model of Hsu et al., 21 one of the diameters in the model is adjusted as a function of temperature to achieve density representability.

Macromolecules
Article for other CG models, following the same expected trend as the length of the chains is increased. The Flory characteristic ratio is defined by C N = ⟨R 2 ⟩/N C l 2 , where l is the atomistic backbone length and N C is the number of carbon−carbon bonds in the backbone. We calculate C N taking the atomistic value l = 1.53 Å and N C ≈ 2N mono ; the results are shown in Figure 7 and compared with other values reported in the literature. The data taken from Milano et al. 10 correspond to simulations at T = 500 K, the same temperature as our simulations, whereas the CG model of Harmandaris et al. 16,17 was studied at T = 463 K and the CG model of Hsu et al. 21 at T = 450 K. The experimental values for the infinite chain are C ∞ = 8.3, 8.5, and 8.6, respectively. 73,74 In Figure 7, a dashed line indicates the expected value of C ∞ = 8.3 at T = 500 K. As can be observed, the characteristic ratio for our model reaches an asymptotic value ∼8 for systems with N mono > 200.
The flexibility of our model is assessed by calculating the angular distribution of the angles formed by three consecutive backbone segments and by backbone−backbone−branch triads as well as for the dihedral angle, defined by as well as for the dihedral angle, defined by branch−backbone−backbone− branch tetrads. The systems with the shortest and the largest chains (N mono = 10 and N mono = 300) are tested, resulting in virtually the same distributions. The results for the N mono = 300 system simulated at T = 500 K are shown in Figure 8. The angular distribution for the backbone ranges from 60°to 180°, presenting two peaks: one at 67°and a second at 115°. These values are considerably smaller than those for CG models with similar structure and connectivity, reported in refs 21, 22, and 25, where a bimodal distribution with peaks at ∼125°and 170°w as obtained from atomistic simulations and used to adjust the bonded interactions of the corresponding CG models. A similar trend is observed for the bbn−bbn−bch angular distribution. The dihedral distribution for our model presents all possible angles with the same probability, while a peak at ∼105°is reported by Hsu et al. for their atactic PS CG model. Our topdown model does not include angular restrictions, which explains the small angles explored by the backbone and the lack of structure in the dihedral distribution. We emphasize again that our model is not developed to resolve the tacticity of PS molecules and is appropriate for the representation of atactic PS. The end-to-end distances and radii of gyration, described earlier in the section, are nevertheless similar to other CG models in the literature. The agreement of the structural properties at one scale but not at the other could be related to the chain flexibility but also to the bulky representation of the backbone segments.
To evaluate the dynamic behavior of the PS model, the selfdiffusion coefficients for the centers of mass (CM) are determined as using simulation trajectories of the melt systems described in Table 3. The CM self-diffusion coefficients as a function of the Figure 13. Snapshots of PS2030 + n-heptane at several temperatures and compositions. The thermodynamic states are color-coded and correspond to those marked with small circles in the temperature−volume fraction phase diagram (cf. Figure 12). The solvent molecules are represented as a continuous transparent material to facilitate the visualization.

Macromolecules
Article chain length N mono are shown in Figure 9. The values of the diffusion coefficient are of the order of 10 −6 cm 2 /s for the N mono = 10 system and of 10 −8 cm 2 /s for the N mono = 300 system. These values are between 2 and 3 orders of magnitude higher than the experimental values, according to refs 75−78. This factor in the time scaling is typical in CG models, where the reduction of friction due to the coarse representation allows faster dynamics. Regarding the scaling with N mono , we observe a single trend in the diffusion. Correlating the data to a powerlaw function, D ∼ N mono a , gives an exponent a = 1.23 ± 0.02, which is closer to the Rouse model (where a = 1) than to reptation behavior, which is characterized by an exponent of a ∼ 2, which is expected for the long-chain systems. This type of two-regime behavior has been exhibited with other models, for example, those in refs 10, 15, and 17. Further understanding of dynamic properties is beyond the scope of our current work; a detailed study of these properties for a CG polystyrene model is discussed by Harmandaris and Kremer in ref 17.
To study the effect of the temperature on the melt density, systems with N mono = 10 are simulated under NPT conditions at P = 1 bar. In Figure 10, we present the resulting densities as a function of temperature including data for other models reported in the literature. We note that the CG model of Hsu et al. 21 features a variable diameter, which was adjusted to be a function of temperature in order to reproduce the experimentally observed density. Problems relating to the representability of the CG models are evident with the CG model of Carbone et al.; 27 the temperature dependence of the density follows an incorrect trend, matching the full-atom density (also shown) only at T = 500 K, the temperature used to develop the model. 27 Our model leads to densities of the order of the expected values; although the trend with temperature is correct, the slope of our data is seen to be larger than that of the other models. As a result, the thermal expansion coefficient α = V −1 (∂V/∂T) P for our system is α = 10.4 × 10 −4 K −1 at T = 300 K, while Hsu et al. report α = 5.50 × 10 −4 K −1 and α = 5.94 × 10 −4 K −1 for the experimental and fully atomistic simulation values, respectively, at the same temperature.
We perform further simulations under NPT conditions to test the effect of pressure, where now the temperature is fixed at a value of T = 500 K. The resulting densities are shown in Figure 11 as a function of pressure, together with CG and atomistic simulation results of Carbone et al. 27 at the same conditions. There are again issues of representability with the model by Carbone et al.; the full-atom density is reproduced only at the pressure of P = 1 atm, which was used to develop the model. 27 Our model gives densities of the order of the expected values and the density−pressure curve has a small slope indicating a high incompressibility. The isothermal compressibility K T = V −1 (∂V/∂P) T for our system is K T = 1.64 × 10 −6 kPa −1 , while Carbone et al. report   Figure 12). The snapshots of the simulation box represent equilibrium configurations of the system at T = 310, 420, and 535 K. (b) Normalized frequencies of the n sol /n tot ratio (cf. Figure 14). The system composition is PS mass % = 24.8 (Φ = 0.18), and the histograms are calculated using d cell 3 = 216 subcells.

Macromolecules
Article K T = 9.3 × 10 −7 kPa −1 and K T = 1.6 × 10 −5 kPa −1 as the fully atomistic and CG simulation values, respectively, at T = 500 K. The description of other quantities, such as glass transition and tacticity effects, are beyond the scope and the objectives of our current work.
Mixtures. The thermodynamic conditions of the polystyrene solutions studied here by molecular-dynamics simulations are summarized in Table 4. From these simulations we present first the results for PS + n-heptane. In Figure 12, we provide the temperature−volume fraction diagram for PS2030 + n-heptane, incorporating the experimental cloud curves given in Figure 1b of ref 39; the bounded regions of liquid−liquid demixing are indicated in the diagram by the gray shading. Volume fractions of the system are defined in relation to the relative values of the segment diameters σ; details are given in the Appendix, exemplified for the case of a PS4800 + n-heptane system. We note that alternative definitions based on the liquid-phase density have little impact on the results. This system is used to adjust the PS−heptane unlike interactions, through simulations at a composition of 24.3 mass % PS (corresponding to a volume fraction of Φ = 0.18).
Simulations are carried out to examine states across a broad region of the phase diagram, indicated by the circles in Figure  12. Snapshots of equilibrium configurations that are representative of these 18 states are shown in Figure 13. All of the snapshots presented in our work are produced using VMD. 79,80 We note that particular care is required in equilibrating the simulations corresponding to the higher-temperature states since the proximity of the critical point of n-heptane causes pronounced fluctuations in the system. Demixing can be easily observed in the snapshots for the extreme cases in temperature, particularly in the lowest-temperature state T = 275 K, for all compositions. However, for intermediate temperatures, it is difficult to distinguish the condition of the mixture. To quantify the state (mixed/demixed) of the mixtures, we follow the approach of Gelb and Muller. 81 Histograms of the local composition of the simulation box are collected. Initially, the box is divided into d cell 3 subcells, where d cell is adjusted for each system, depending on the composition, to give the best possible noise-to-signal result. The total number n tot of segments in each subcell as well as the segments n sol corresponding only to solvent are accumulated in histograms of the n sol /n tot ratio, fixing the width at 0.02 for all cases. This is repeated over hundreds of equilibrated configurations to provide good statistics. Finally, the frequency f for each n sol /n tot is normalized with respect to the number of configurations and subcells, so that ∑ i f i = 1. A symmetrical distribution centered around the global composition of N sol /N tot indicates that, on average, the local compositions are similar throughout the system. By contrast, demixing is characterized by a bimodal distribution as its fingerprint, indicating how some subcells have a predominantly solvent-rich content, while others are mainly populated by PS segments.
The normalized distributions for each thermodynamic state are presented in Figures 14a−c. The vertical line indicates the overall system composition ratio N N / sol tot (in terms of number of segments), for each simulation box. From Figure 14, at each composition, we can observe that the system at the lowest temperature of T = 275 K is characterized by a bimodal, highly asymmetrical distribution, indicating that T = 275 K corresponds to a demixed state. This is satisfying, given that T = 275 K is below the reported UCST (311 K) 39 for this system. At higher temperatures, the distribution gradually transforms into a symmetrical bell shape; this is particularly apparent for the systems at T = 400, 450, and 500 K. This behavior is in agreement with the expectation from the fluidphase diagram presented in Figure 12. At the highest temperature of T = 535 K the asymmetry appears again at each composition; this corresponds to a temperature above the reported LCST (515 K) 39 for this system. The histograms at T = 325 K exhibit an intermediate behavior, which can be attributed to the proximity to the limiting UCST (see Figure  12). Exploratory evaluations of the Hsu et al. 21 model showed that it predicts the UCST at higher temperatures than the SAFT-γ model, inconsistent with the experimentally observed phase behavior. This occurs despite slightly better predictions of the temperature−density behavior with the Hsu et al. 21 model, suggesting an advantage of the top-down approach advocated here to describe phase separation.
We test the representability of the unlike-interaction parameters by exploring the behavior of solutions comprising the longer polymer, PS4800. The results at the thermodynamic conditions indicated in Figure 15a correspond to the histogram distributions presented in Figure 15b. In this case, the LCST is situated at a lower temperature of T = 477 K, 39 so that our highest-temperature simulation is sufficiently above the LCST to capture the demixing expected in the high-temperature region. The histogram distributions for these systems are asymmetrical for the lowest and highest temperature, T = 310 K and T = 530 K. On the other hand, the distribution is very symmetrical for the case of T = 420 K, while the distributions for T = 370 K and T = 465 K exhibit an intermediate behavior, in agreement with the proximity of these Figure 16. (a) Temperature−volume fraction phase diagram for PS2030 + n-hexane (cf. Figure 12). The snapshots of the simulation box represent equilibrium configurations of the system at T = 280, 400, and 495 K. (b) Normalized frequencies of the n sol /n tot ratio (cf. Figure 14). The system composition is PS mass % = 24.1 (Φ = 0.18), and the histograms are calculated using d cell 3 = 216 subcells.

Macromolecules
Article thermodynamic states to the experimental LCST and UCST reported by Cowie and McEwen. 39 The phase behavior for PS + n-hexane is more challenging to reproduce since an "hourglass" phase diagram is expected in one of the systems but not in the other, as can be seen in Figures 16a and 17a. The unlike-interaction parameters in this case are adjusted with simulations of the system PS2030 + nhexane at a composition of mass % = 24.1, as described in our earlier subsection Unlike Interactions for Mixtures. As for PS + n-heptane, simulations are carried out to examine the system across a broad region of the phase diagram; the thermodynamic states examined are indicated by the circles in Figure 16a. The resulting composition distributions are shown in Figure 16b. Again, the lowest-and highest-temperature states are characterized by asymmetrical bimodal distributions, indicating a demixed system. This can be appreciated from the snapshot of the configuration at the lowest temperature studied, T = 310 K, which corresponds to a state located far from the UCST. For this system, the practical upper limit of temperature for simulations is T = 495 K. Simulations at higher temperatures are not meaningful due to the fluctuations caused by the proximity to the critical temperature of n-hexane. The system at T = 330 K yields a histogram with similar asymmetry to the highest-temperature system. This behavior is explained by the proximity of this state to the limiting UCST.
As a second assessment of the PS−hexane unlike-interaction parameters, the same force fields are used to simulate the PS4800 model. The experimental "hourglass" phase diagram, 39 and the thermodynamic states selected for simulation are depicted in Figure 17a, including snapshots corresponding to the relevant equilibrium configurations. The local-composition histograms are shown in Figure 17b,c. The systems at the lower PS composition, mass % = 24.6 (Φ = 0.18), exhibit asymmetrical distributions at all of the simulated temperatures, as illustrated in Figure 17b. The histograms for the higher PS composition, mass % ∼ 49.7 (Φ = 0.40), are represented in Figure 17c. The states at T = 280 K and T = 490 K are characterized by very asymmetrical bimodal distributions, with strong peaks at low n sol /n tol and at n sol /n tol = 1, indicative of liquid−liquid demixing. The state at T = 400 K, on the other hand, is characterized by a distribution that although slightly asymmetrical, does not feature a second peak; this distribution resembles the bell shape indicative of complete miscibility corresponding to a single-phase system. It is particularly gratifying that, collectively, these results are entirely consistent with the "hourglass" phase diagram, illustrated in Figure 17a, observed experimentally by Cowie and McEwen. 39

■ CONCLUSIONS
We present a coarse-grained model for polystyrene and its solutions in aliphatic solvents, exemplified by n-hexane and n-heptane. The force-field parameters are obtained using the SAFT-γ Mie EoS. 32,43 This top-down approach is a simple way to generate CG models of long or complex molecules.
Systems of pure solvents are studied by MD simulations and theory, with good agreement between the phase diagrams generated using both methodologies and those from experiment. Melts of the PS models are tested under NPT and NVT conditions to study the effect of the polymer chain length and of temperature and pressure. The quality of the results, in terms of the density, end-to-end distance, and radius of gyration, is similar to that obtained using other CG models reported in the

Macromolecules
Article literature for atactic polystyrene, generated via bottom-up procedures. Our CG polystyrene is simple enough to be used in simulations of large systems, such as those required to study PS in solution. At the same time, the description maintains the main features needed to provide good representability, i.e., a reliable PS model that can be used under different thermodynamic conditions. In the Introduction we alluded to the general issues of transferability (the ability to transfer a moiety (segment or bead) corresponding to a particular functional group in one molecule to another molecule containing the same functional group), representability (the ability of a model to represent the parent experimental system at different thermodynamic states), and robustness (the ability of a model to provide reliable descriptions across a wide range of different properties). We have demonstrated transferability of our CG PS model through the study of PS at significantly different molecular weights. We have also demonstrated representability of our models, through simulations both of pure-component and mixtures over broad regions of their respective phase diagrams. In this context we note, however, that the transferability of our mixture models to systems which are markedly different to those we have considered, and the representability for conditions differing very widely from those considered, have not been thoroughly tested. Any generalization of the assignment of the k ij values to mixtures or conditions dramatically different to those considered in our current work should be performed with some care. One would expect the results to be very sensitive to the fine details of the values of the binary-interaction parameters. The fluid-phase behavior showcased here is the result of a very fine balance between entropic and energetic contributions.
We have demonstrated the robustness of our models in terms of structural and thermodynamic properties. Since our models are not informed by atomistic simulations, they lack some inherent structural information. This, however, is not a limitation of the top-down approach. SAFT models can be "decorated" with appropriate intramolecular interactions without compromising their accuracy in describing macroscopic properties. 48 This is particularly important in modeling transport properties such as diffusion and viscosity. An investigation of this in the context of our PS model was beyond the scope of our current study and will form the basis of future work.
A natural extension of the present work in relation to mixtures involving PS would be the exploration of more complex systems. Simulation of multicomponent mixtures, over broader ranges of thermodynamic conditions, is an attractive possibility that would improve our understanding of the interplay of effects between the thermodynamic variables and the properties of the system for cases where the information is not accessible to experiment or theory.
where N chain is the number of PS chains, N mono is the number of styrene units, and σ is the backbone (bbn), branch (bch), or solvent (sol) diameter. This definition is exact for hard (athermal) systems, but not for soft spheres. We recognize that a more general definition of the volume fraction One should note that this last definition relies on the knowledge of the temperature-dependent densities and hence is impractical. Numerically, however, both definitions produce similar results. As an example, the PS4800 + n-heptane system has a fixed volume fraction of Φ = 0.18 following eq 7. To calculate the corresponding values through eq 9, we use the temperature-dependent PS densities obtained with the pure melt simulations, N mono = 10 system, and the densities of pure n-heptane. The resulting values range from 0.17 to 0.16 for temperatures between 310 and 465 K, while only the highesttemperature system (at 530 K) has a volume fraction that deviates from the fixed value by more than 25%, with a volume fraction of Φ = 0.13.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.macromol.6b02072. Simulation details of solvents, PS melts, and PS−solvent mixtures; plot of the PS−PS contribution to the system energy as a function of time (PDF)