Using the SAFT-γ -Mie to Generate Coarse-Grained Force Fields Useable in Molecular Dynamics Simulations: Describing the Micellar Phases of Polyalkylglycols in Aqueous Solutions

: In this work, we showcase the potential of the SAFT-γ -Mie equation of state to obtain, in a top-down approach, coarse-grained molecular models capable of describing micellization phenomena. For this purpose, we selected polyethylene oxide (PEO)/polypropylene oxide (PPO) triblock copolymers in aqueous solutions as a case study. The model was fitted to reproduce the thermophysical properties � liquid densities and vapor pressures � of small PEO/PPO oligomers as well as the activity coefficients or liquid − liquid phase diagrams of their binary mixtures with water. At low concentrations, the MD simulations show the formation of spherical micelles of morphology comparable to experiments, with the characteristic star-shaped or flower-shaped morphology depending on the polymer block sequence. The local structure and molecular conformations were analyzed, showing distinctions between the two types of aggregates. At higher concentrations, depending on the co-polymer sequence, the molecules form larger micelles or coalesce into larger aggregates when signs of phase separation are observed, as expected.


INTRODUCTION
−4 At its core, the SAFT EoS describes the interactions between molecules with an explicit potential, which enables a higher reliability than other equations of state.The use of a group contribution EoS, such as SAFT-γ-Mie, 5 conveys an easy way to obtain transferable coarse-grained force fields for molecular simulations based on an accurate description of thermophysical properties.Even though the SAFT EoS can be used to calculate thermodynamic equilibrium properties, it lacks, on its own, the ability to predict dynamic properties or structural ones, such as self-assembly phenomena.However, in combination with molecular dynamics (MD) simulations, it is possible to study the supramolecular organization of complex systems that undergo phase transitions, such as liquid crystals 6,7 or micellar solutions, 8 and also study interfacial 9−11 or transport properties. 12This transition of the SAFT model into molecular simulations has its caveats, though, as referenced by Muller and Jackson. 13The inclusion of associative interactions is not trivial given the discontinuous nature of its treatment, and the shape-factor, which represents the weight of a segment in the calculation of the properties with the SAFT EoS, needs to be a whole number, thus restricting the parameter phase space.The SAFT-γ-Mie force field for simulation has been used with great success to predict the solubility of various apolar polymers in hydrophobic media. 14,15However, little work has been done regarding aqueous solutions of polymers, possibly due to the difficulty of capturing the unique properties of water, especially when the association term is absent.Nonetheless, preliminary work on the self-assembly of C i E j surfactants in aqueous solutions was disclosed by Muller and Jackson. 13n this work, we attempt to apply the SAFT-γ-Mie EoS in a top-down approach, to model the micellar phases of a family of triblock polymeric compounds containing polyethylene oxide (PEO) and polypropylene oxide (PPO) monomers, known as poloxamers.−25 Poloxamers behave as amphiphiles due to the difference in hydrophobicity of the PEO and PPO moieties, and, depending on how the blocks of the polymer are arranged, the properties of the selfassembled entities can differ significantly.For PEO-PPO-PEO sequences, commonly known as L-type poloxamers, the structure of the micelles is star-shaped.On the other hand, PPO-PEO-PPO sequences are known as R-type poloxamers, and their micelles are flower-shaped.In both cases, the hydrophobic PPO groups populate the core of the micelles, while the hydrophilic PEO segments occupy the outer shell.Herein, we showcase the potential of the SAFT-γ-Mie EoS to be used in a top-down approach to model the aggregation of complex polymeric systems in water solutions.A group contribution methodology was used to capture the thermodynamic properties of short oligomers and then applied to describe the aggregation properties of polymeric surfactants in aqueous solutions using molecular dynamics simulations.A new coarse-grained force field, including the intermolecular and intramolecular potentials, for the description of poloxamers is proposed here.

SAFT-γ-Mie EoS.
The SAFT is one of the most successful molecular-based equations of state used for the description of nonideal systems, such as associating molecules.SAFT is usually written in terms of the Helmholtz free energy, A res , where different terms are added to take into account the interactions of the molecules integrating the fluid.For a system of associating chains, the residual Helmholtz free energy can be described by where A monomer , A chain , and A assoc stands for the free energies of the reference (monomer or segment) fluid, the formation of the chains, and the association contributions, respectively.In order to obtain a model suitable for MD simulation, the association term was not included in this work.In the case of SAFT-γ-Mie, the molecules are represented by tangentially linked spherical monomers that interact via the Mie potential.
For two particles i and j separated by a distance r, the potential can be written as follows where ε ij is the potential well depth, σ ij is the characteristic distance where the potential is null, and λ r /λ a are the repulsive and attractive exponents.In this work, the attractive exponent was set to 6, following London's definition for dispersive interactions.C is a prefactor that is solely dependent on the exponents, and can be written as follows = i k j j j j j j y { z z z z z z i k j j j j j j y { z z z z z z The cross-interaction parameters between different bead types i and j can be calculated from the same type interactions using the Lorentz−Berthelot combining rules written below.
In this work, σ ij and λ r,ij were not fitted, while in the case of ε ij , the cross-interaction was tuned via the adjustable parameter, k ij , to allow deviations from the geometric mean.All the SAFT-γ-Mie calculations were made using the python module SGTPy, made available by Chaparro Maldonado. 26.2.SAFT-γ-Mie Force Field.To fully define the force field, in addition to the intermolecular interactions underlined in the previous section, the intramolecular bonded interactions also need to be specified in order to realistically capture the internal degrees of freedom of the polymer chains.An atomistic force field (OPLS-AA 27 ) was used to reference the intramolecular configurations in order to obtain the harmonic bond and angle potentials, U bond and U angle , for the coarsegrained model.
In eqs 7 and 8, k b and k θ are the harmonic constants, while r 0 and θ 0 are the equilibrium distance and angle, respectively.

SAFT-γ-Mie Building Blocks.
Poloxamers are triblock co-polymers consisting of PEO and PPO segments linked together, and their structural formula is represented in Figure 1.In this work, two poloxamers of similar molecular weight (∼2000 Da), L-type L35 and R-type 10R5, were selected.
The scaling from atomistic to coarse-grained grouping followed a 3:1 and 4:1 mapping to represent the geometry of the PEO and PPO groups, respectively.This ratio provides a good balance in order to obtain a model that captures the interactive nature of the polymer molecules, while also allowing for a considerable speed-up in run time compared to atomistic models.Furthermore, it also allows for good transferability and scaling due to the fact that each bead of the model corresponds exactly to one monomeric unit of the polymer.As illustrated in Figure 2, three different bead types were developed to model the polymer compounds: an end group corresponding to the CH 3 OCH 2 segment, one PEO group that represents the CH 2 OCH 2 monomer, and a PPO bead that corresponds to the CH(CH 3 )OCH 2 or CH(CH 3 )-OCH 3 unit.The fitting procedure is described in section 3.1 in more detail.
The water model used in this work was developed by Lobanova et al. 28 and incorporates two water molecules per bead.It was developed to reproduce the liquid density and surface properties of water over a large range of temperatures.A crucial advantage of this model is that the melting temperature is closer to the experimental value than the other models described in the work of Lobanova et al.

MD Simulation Details.
Molecular dynamics simulations were carried out for 200 polymer molecules, L35 or 10R5, in cubic simulation boxes containing the corresponding water molecules to satisfy total polymer concentrations of 1 and 5% wt.The initial configurations containing randomly inserted molecules were generated with packmol. 29The initial configuration underwent an energy minimization stage followed by a compression/equilibration run before the production run.Simulations were made with Gromacs 2018.1 software 30 using the leap-frog algorithm to integrate the equations of motion.A timestep of 10 fs was used, and the cut-off distance was set to 1.4 nm, which corresponds to three times the largest value of σ.Beyond the cut-off distance, dispersion interactions were corrected by adding the missing energy and pressure to the system.All simulations ran for a total of 1000 ns in the NpT ensemble.The pressure and temperature were kept constant with the Parrinello-Rahman 31 and Nose-Hoover 32,33 barostat and thermostat, with coupling constants of 10 and 0.5 ps, respectively.Pressure was set to 1 bar, and the temperature was kept at 298.15 or 343.15K.The algorithm for determination of aggregation has been described previously, 34 where the cut-off distance between PPO beads was set to 0.7 nm, the first minimum of the PPO-PPO radial distribution function.

SAFT-γ-Mie Parameterization.
The SAFT EoS was fitted to describe the vapor pressures and liquid densities of small oligomer molecules containing ethylene oxide and propylene oxide groups.The parameters for the pure oligomers   35 were obtained without the need to tune unlike interactions, under the assumption that the cross-interactions between ether groups follow the combination rules due to their similar chemical nature.Diglyme, triglyme, tetraglyme, and dipropylene glycol dimethyl ether (DPGDME) were the chosen compounds for which experimental data was found in the literature.A total of nine parameters were fitted, corresponding to three sets of same-type parameters (ε, σ, and λ r ) for each group.This differs from the approach previously employed by Lobanova 8 to model CH 2 −O−CH 2 groups, in which the repulsive exponent was fixed as λ r = 19.The shape factor, s k , was set to 1, and the number of monomers, m, corresponds to the number of beads of each type within the molecule.The densities of the glyme compounds were previously measured as a function of temperature and pressure, 35 thus providing enough data to describe the properties of these compounds.In the case of DPGDME, the available data is scarcer and was mainly sourced from data sheets provided by manufacturers.Overall, the fitted properties (see Figure 3a,b) reached a very good agreement with the experimental results, with average absolute deviations (see Table 1) below 10% for vapor pressures and below 1% for the densities, which is well within the expected deviations obtained by other authors using similar approaches. 12,15Furthermore, the derivative volumetric properties of the polyether compounds (see Supporting Information) are also very well described.Calculated isothermal compressibilities deviate less than 6% from the experimental data, while the isobaric expansion coefficients have average absolute deviations below 10% for all compounds.With an increase in the chain length of the oligoether compounds, there is no clear trend for higher deviations for the experimental data, thus boding well for when scaling to higher molecular weight polymers is required.
The cross-interaction parameters between the polymer groups and water were obtained by fitting the potential well depth, ε ij , in order to reproduce the water activity coefficients (in the case of the glyme compounds) or the LLE equilibrium diagram (for DPGDME).The water activity coefficients were chosen instead of the VLE data due to the inability of the chosen nonassociative water model to accurately describe the saturation pressures. 28The results in graphical form can be seen in Figure 4. Previous SAFT approaches using an association scheme and with fewer restrictions in the parameter phase space, despite the excellent agreement of the VLE curves with the experimental data, did not perfectly capture the behavior of the water activity coefficients. 38onetheless, it is reassuring that the results obtained in this work for the water activity coefficients follow the same trend obtained by Crespo et al., 38 where an increase in the saturation pressure leads to lower deviations from ideality.In the case of the LLE diagram calculated with the SAFT EoS, the lower critical solution temperature (LCST) behavior could not be captured, despite some qualitative agreement in the two-phase zone above the LCST.Attempts to fit the softness of the Mie potential by changing the repulsive exponent were unable to improve the description of the phase diagram, and as such, only the well depth cross-parameter was tuned.Tables 2 and 3 contain the full set of parameters obtained for the SAFT-γ-Mie EoS.

Intramolecular Potentials.
The description of the intramolecular interactions can have a significant impact on the structural and dynamic properties of molecular systems, as previously demonstrated in the case of SAFT-γ-Mie derived force fields. 40As such, following a similar methodology to the one proposed in the aforementioned article, we used atomistic simulations as a reference for the internal configurations of the polymer molecules.The OPLS-AA force field was chosen to model oligomers comprising a total of 8 monomers of either PEO or PPO groups.Details on the atomistic simulations used for the intramolecular potential parametrization can be found in the Supporting Information.
Atoms corresponding to the same CG bead were grouped, and the distance/angle distributions were calculated between the centers of mass of linked groups.The harmonic potentials (eqs 7 and 8) were fitted to reproduce the potential energy obtained from the MD trajectory after applying the Boltzmann inversion equation where x is either the bond length, r, or angle, θ, between linked CG beads, RT is the ideal gas constant multiplied by the temperature, P x is the distribution probability, and K is an arbitrary constant that can be used to shift the energy profile.
The calculated potentials can be seen in Figure 5, and the set of fitted parameters can be found in Table 4.It should be noted that the equilibrium bond length obtained from the fitting does not correspond to the ones used in the simulation force field due to the definition of the SAFT EoS according to tangentially bonded segments.Thus, the value of r 0 included in the definition of the force field is equal to the value of σ.In general, the distribution of distances and angles between PEO groups shows a single peak and is well represented by the harmonic potential.As for PPO groups, due to the ramified nature of the polymer, rotation across the bonds of the main chain results in bond and angle distributions of higher complexity, and thus, the agreement with the harmonic potential fit is not as good.Nonetheless, the key characteristic differences between the two types of segments are well represented: the higher rigidity and length of the PPO−PPO bonds and the lower average angle of the PPO−PPO−PPO segments.

MD Simulations.
The new SAFT-based force field developed in this work was used to perform CG-MD simulations of aqueous solutions of L35 and 10R5 at 1 and 5% wt concentrations and 333.15 or 298.15K temperatures.Simulation snapshots after 1000 ns of simulation time are displayed in Figure 6.At the lowest concentration, L35 and 10R5 formed the archetypical star and flower-like micelles, respectively, also found in the work of Peŕez-Sańchez 23 using the MARTINI 2 force field.The inner and outer shells are segregated, as the inner core is composed of PPO segments, and the outer shell is occupied by the PEO moieties, following the different affinities towards water.As the concentration of L35 increases from 1 to 5% wt, attractive interactions between the PEO groups of L35 neighbor micelles promote micelle growth (see Figure 6a,b) by means of micelle coalescence and Industrial & Engineering Chemistry Research monomer absorption.At 1% wt, both L35 and 10R5 exhibited similar micellar distributions, but as soon as the concentration of 10R5 is increased to 5% wt, an increased aggregation process was observed, hinting towards signs of phase separation.This different behavior between L35 and 10R5, with an increase in concentration, is supported by the cloud point divergence observed experimentally, 23 given that at a temperature of 333.15K and a concentration of 5% wt, we are below the LCST for L35 but above the LCST for 10R5.Simulations also show that, at higher concentrations, intermicellar contacts are seen between L35 PEO segments, with some interdigitation between neighboring aggregates.Conversely to previous simulation studies, 23 physical crosslinking is not prevalent for any of the systems studied.
With an increase in concentration from 1 to 5%, the average aggregation numbers of L35 increase from 24 to 32 with micelles homogeneously distributed along the simulation box.On the other hand, 10R5 goes from aggregates containing an average of 20 molecules at 1% wt to full clustering at 5% wt.Intramicelle distance distributions in L35 and 10R5 solutions at 1% wt are presented in Figure 7, allowing a quantitative evaluation of the structure as well as the size of the aggregates.The distribution of the PPO groups does not differ between micelles comprising L35 or 10R5, meaning that no significant  Industrial & Engineering Chemistry Research structural differences in the inner core of the two types of micelles are observed.However, the wider distribution of PEO groups in L35 indicates that L-type poloxamers, where the ends are composed of PEO, tend to be more readily hydrated by the water subphase than R-type counterparts.This can be easily visualized in the simulation snapshot insets of Figure 7, where 10R5 chains tend to be more crumpled, while L35 segments are mostly extended.This observation about the way the PEO groups are hydrated helps to rationalize why, at a higher concentration of 5 wt %, only the 10R5 polymer undergoes phase separation.
At a 5% wt concentration, a decrease in temperature from 333.15 to 298.15 K does not result in significant differences.In the case of L35, spherical micelles are still observed with comparable aggregation numbers.As for 10R5, a clear segregation between water and the organic phase is present, but with some interesting features.At a higher temperature, 333.15 K, due to the thermal agitation, the prolate-shaped aggregate is more disorganized, and the number of dispersed 10R5 molecules in the aqueous phase is higher than at a lower temperature of 298.15 K.
Radial distribution functions (RDF) can provide detailed insight into the structure of the L35 and 10R5 solutions.Figure  8a,b show the RDF between alike polymer chains for the systems at 333.15 K and 5% wt concentration.Despite the significant difference in the aggregation behavior of the two polymers, the short-distance interactions are very similar, given the identical amplitude of the first and second order of coordination peaks.The major differences between the two radial distribution functions are seen at long-range distances due to the discrepancy in size of the respective aggregates.Despite the similarity between the RDFs of the polymers when using the whole molecule as a reference, there are some significant differences if we consider the PEO/PPO groups separately.These radial distribution functions are presented in  Figure 8c for all PEO/PPO combinations.Given the lower probability of PEO-PPO contacts, we can infer that the segregation between PEO and PPO groups is higher in the system containing L35 molecules.This directly relates to the higher organization of L35 molecules within the spherical micelles, where PPO and PEO regions are well established, whereas 10R5 molecules are not as ordered in the immiscible organic phase.
The end-to-end distance distribution analysis for the two polymers at 1% wt concentration is exhibited in Figure 9.It is interesting to see that, despite the similarity between the chemical composition, aggregation numbers, and sizes of the micelles of the two polymers, the end-to-end distributions are rather different.L35 molecules, which contain PEO ends, are, on average, much more extended than the reverse 10R5 polymer, whose distribution is much narrower.Even though the two molecules contain the same segments, the order in which they are bonded heavily influences the internal conformations of the molecules and, thus, their potential energy.This difference in the conformations of the two polymers is illustrated in Figure 9, where it is clear that the PPO−PEO−PPO configuration in 10R5 moieties tends to form loops and, thus, flower-like micelles, while the star-like micelles contain mostly extended molecules in star-like assemblies.The bends on the 10R5 molecules contribute to the lower stability of the micellar phase when compared to the L35 counterpart, which might partially explain the higher tendency of 10R5 to phase separate.

CONCLUSIONS
In this work, a coarse-grained force field for polyalkylglycol molecules was developed in a top-down approach using the SAFT-γ-Mie EoS to obtain the intermolecular interaction parameters.The fitting process, applied to small oligomers containing PEO and PPO groups, produced good quantitative results for the pure compounds, in good agreement with experimental results.The cross-interaction parameters were tuned to reproduce the polymer-water binary mixture properties, namely the water activity coefficients and liquid−liquid diagrams.The activity coefficients followed the experimental data with good accuracy, while the calculated LLE diagram of DPGDME + water was unable to capture the existence of a lower critical solution temperature.
The supramolecular organization of triblock co-polymers of PPO and PEO was tested with the proposed force field for polymers of different block sequences using molecular dynamics simulations.At a concentration of 1% wt in polymer, the model is able to predict the micellar phases of the poloxamers in water media.At this concentration, the two types of polymers do not differ in size but do differ in morphology, as L35 adopts a star-like shape and 10R5 a flower-like shape.At a higher concentration of 5% wt, the aggregation behavior of the two polymers diverges, as the 10R5 phase separates while L35 maintains its micellar phase with larger aggregates.The conformational stability is appointed as a possible reason for the different aggregation behavior of the two polymers.A decrease in temperature does not have a significant impact on the aggregation of L35; however, the 10R5 assemblies, due to lower thermal agitation, form more cohesive aggregates.

Figure 1 .
Figure 1.L-type and R-type poloxamer triblock co-polymers general formulas are shown on the left and on the right, respectively.For L35, x = 11 and y = 16, while for 10R5, x = 22 and y = 8.

Figure 2 .
Figure 2. Coarse-grained (CG) mapping considered for PEO (left) and PPO (right) segments.The colors represent the three different types of beads used in the model.

Figure 3 .
Figure 3. Vapor pressures (a) and liquid densities of diglyme (b), triglyme (c), and tetraglyme (d) fitted with the SAFT-γ-Mie equation (lines) compared to the experimental results (symbols).Vapor pressures were sourced from refs 36 and 37. Densities were taken from the work of Navarro. 35

Figure 4 .
Figure 4. Binary equilibrium properties for oligomer + water systems, where experimental data are displayed as markers andSAFT-g-Mie calculations are shown as solid lines.Water activity coefficients in systems containing diglyme (a), triglyme (b), and tetraglyme (c), plus liquid− liquid equilibrium for the DPGDME + water system.Experimental data can be found in the works of Crespo38 and Zhao.39

Figure 5 .
Figure 5.Bond stretching (a) and angle bending (b) potentials fitted to reproduce the distribution obtained from atomistic simulations using the OPLS-AA force field.

Figure 7 .
Figure 7. Probability distribution of the PPO and PEO segments to the center of mass (CoM) of the largest observed micelle in the simulation at 333.15 K. (a) L35 at 1% wt; (b) 10R5 at 1% wt PEO and PPO beads are colored pink and green, respectively.

Figure 8 .
Figure 8. Radial distribution functions between polymer molecules (a), polymer-water (b), and combinations of PEO/PPO groups (c) at 333.15 K and concentration of 5% wt.

Table 1 .
Absolute Average Deviations (AAD) for the Calculated Properties upon Fitting of the SAFT-γ-Mie EoS

Table 3 .
Fitted Cross-Interaction Parameters a a All cross-parameters not referenced on the table below were calculated using the combination rules (eqs 4−6)