Multiscale Modeling of Binary Polymer Mixtures: Scale Bridging in the Athermal and Thermal Regime

Multiscale Modeling of Binary Polymer Mixtures: Scale Bridging in the Athermal and Thermal Regime

J. McCarty    M.G. Guenza111Author to whom correspondence should be addressed. Electronic mail: Department of Chemistry and Institute of Theoretical Science, University of Oregon, Eugene, Oregon 97403
July 3, 2019

Obtaining a rigorous and reliable method for linking computer simulations of polymer blends and composites at different length scales of interest is a highly desirable goal in soft matter physics. In this paper a multiscale modeling procedure is presented for the efficient calculation of the static structural properties of binary homopolymer blends. The procedure combines computer simulations of polymer chains on two different length scales, using a united atom representation for the finer structure and a highly coarse-grained approach on the meso-scale, where chains are represented as soft colloidal particles interacting through an effective potential. A method for combining the structural information by inverse mapping is discussed, allowing for the efficient calculation of partial correlation functions, which are compared with results from full united atom simulations. The structure of several polymer mixtures is obtained in an efficient manner for several mixtures in the homogeneous region of the phase diagram. The method is then extended to incorporate thermal fluctuations through an effective parameter. Since the approach is analytical, it is fully transferable to numerous systems.

I Introduction

Mixtures of polymers with different compositions are of great interest for many technological and industrial applications. For example, multicomponent polymer mixtures facilitate the custom tailoring of desired physical properties, leading to the creation of new materials with enhanced performance.Sperling () However, since new materials are typically processed in their liquid phase, the ability to predict the mixing behavior of various polymer composites is a desirable computational task. Because of the underlying phase diagram, polymer mixtures can phase separate at specific thermodynamic conditions of temperature and polymer volume fraction. The structure and dynamics of each component strongly depends on the region of the phase diagram that is explored, and different length scales characterize the thermodynamic and dynamical behavior of the mixture, depending on how far the system is from its phase separation.

While computer simulations are useful in elucidating the structural and dynamic complexity of these systems, they are limited in the extent of information that can be collected by the precision of the calculation, which degrades with the number of simulation steps, and by the power of the available computational machines.simul () Already for polymer liquids it is not possible to collect all the relevant information from a single simulation when long polymers are involved, as information spreads over many order of magnitude in time and length, going from the monomer length scale, , to the overall polymer dimension of the radius-of-gyration, . The latter scales as the square root of the degree of polymerization, , which can be fairly large for polymers of interest.binderbook () The added complexity of dealing with polymer mixtures arises from the presence of a new length scale of concentration fluctuations, , which diverges as the system approaches its spinodal curve. The big challenge for such systems becomes simulating the mixtures not only at different temperatures and compositions, but also with a box size that increases with the scale of concentration fluctuations. Furthermore, the equilibration step in molecular dynamics (MD) simulations becomes extremely slow, as the time scale in which polymer mixtures phase separate diverges as the temperature approaches the critical point. Finally the timescale of relaxation of the mixture depends on the proximity of the system to the glass transition of its components. In conclusion, classical MD simulations of polymer mixtures have been limited so far to high temperatures close to the athermal limit. In that region the thermodynamic behavior is mainly guided by entropy and local packingbinderbook (), while enthalpic effects may be accounted for within a high temperature perturbative framework.JARAMILLO (); HEINE ()

In order to extend the ability of MD simulations to make quantitative predictions about polymer blends at temperatures approaching the spinodal, information is needed on multiple length-scales of interest. As the system approaches its phase transition and the largest length scale of interest diverges, one might be tempted to neglect the details of the local scale; however, such a procedure risks losing pertinent information, since the monomeric structure of a blend’s component largely determines its glass transition temperature, as well as the shape of its phase diagram. Nonetheless, capturing the global structure is equally as important since an accurate determination of the thermodynamics of the system requires obtaining the total distribution of particles in the ensemble. For example, while system-specific local information is contained within the correlation hole region of the pair distribution function, g(r), the global structure is still required, as truncation of g(r) will prevent accurate Fourier transform to the structure factor where the low k limit is important.ALLEN () Thus, both local and global scale information need to be accounted for to achieve a complete thermodynamic and dynamic description of the mixture.

One strategy to improve our ability to simulate mixtures of complex fluids and collect information on an increasingly large range of length scales is to use a hierarchical approach, such as the multiscale modeling procedure that we present in this paper. In such a multiscaling scheme the system of interest is represented at various levels of coarse-graining and simulations are performed for each of these coarse-grained representations to enhance computational efficiency. In a second step, the information obtained from those simulations is collected and combined into a complete description, bridging information at all the length scales of interest.Plathe ()

For any coarse-grained model to be capable of faithfully reproducing the desired characteristics of the system, the coarse level of description needs to include the essential features of the atomistic structure while averaging out the remaining degrees of freedom. Several coarse-graining procedures have been developed in the literature at various levels of “coarseness”. At the lowest level, the united atom (UA) description, which represents groups of units, with , or into an effective UA along the chain, is a particularly useful representation, for which potentials have been optimized to reproduce the structure for most polyolefins.UApapers () Higher levels of structure-based approaches simplify interactions, while remaining close to the chemical structure, facilitating scale hopping between the different levelsPraprotnik (); however, these methods remain computationally demanding and do not adopt an extreme enough level of coarse-graining to mitigate the limitations of fully atomistic simulations in reaching the large length scales of phase separation.Baeurle () In order to simulate the increasingly large scale phenomena of interest to material scientists, an even coarser level of description is needed, while remaining formally related to the finer structure.

The highest level of coarse-graining in our procedure, represents each polymer chain as a soft colloidal particle, centered at the polymer center-of-mass. Soft colloids interact through a Gaussian-like repulsive potential of the order of the polymer radius-of-gyration, , which defines the characteristic length scale of this description. Large-scale behavior is easily generated in a “mesoscale” simulation of the coarse-grained, mixture and information is collected on a scale equal to or larger than the polymer radius-of-gyration. Recently, we have shown that our mesoscale simulations capture the global structure of the polymer mixture providing information over much of the miscible region of the phase diagram.McCarty1 ()

Much attention has been given recently to relating mesoscopicAkkermans () or field-theoretical model parametersSewell () to their more detailed atomistic representation. For example, Groot and Warren were able to express the repulsive soft potential parameters used in dissipative particle dynamics (DPD) in terms of the Flory-Huggins parameter, providing a bridge between DPD simulations and atomistic simulations, from which solubility parameters can be directly calculated.GROOT () This method has been successfully implemented in mesoscale simulations of biological membranes once the model is parameterized for specific molecular structures.GROOT2001 () In our approach, the effective parameter is included directly in the model as an input parameter which determines the length scale of the of concentration fluctuations. This parameter provides a convenient link between an atomistic or united atom level description and a mesoscale picture, similar to other mesoscale modeling approaches;GROOT (); Sewell (); Titievsky () however, our approach has the advantage over other methods such as DPD in that the soft-core potential depends directly on the parameter and not on any phenomenological expression to bridge the two theories. In other words, once the parameter is specified, our multiscale modeling procedure will produce the expected phase behaviorMcCarty1 (), and the task in practice is to determine this parameter for a given model.

In the current paper, we focus on bridging the length scales of binary homopolymer mixtures of different chemical architecture to provide a complete description of the static structure of the blend under a given set of thermodynamic conditions. Blends of polyolefins represent a good test case for our procedure since united atom parameters have been well optimized for these systems, and they provide a stringent test of the method, as subtle differences in chemical architecture, such as the extent of branching, greatly affects the global structure and miscibility. One advantage of our coarse-graining and multiscale modeling approaches is that since they are based on the formalism of the Ornstein-Zernike equation, the rescaling of the structure to recover the monomer level of description is straightforward. Furthermore, we have shown that the correct dynamics can be obtained directly from simulations of the coarse-grained structures by proper rescaling.Lyubimov () Thus, by performing simulations of mesoscopic particles and subsequently reinserting the relevant chemical details afterwards, we are able to obtain quantitative information about large-scale properties.

Large scale information, obtained form mesoscale MD simulations, is later combined with local-scale united-atom simulations of the same system, using a multiscale formalism obtained from the solution of the Ornstein-Zernike equation, presented in this paper. In a spirit similar to the united atom simulations to which we test our approach, we initially take the high temperature limit as a reference state, and set in order to calculate monomer level total correlation functions. The model is then extended to incorporate thermal effects on the concentration fluctuation through the interaction parameter, , which serves to bridge the gap between mesoscale and atomistic models.

This multiscale modeling approach, is an extension to polymer mixtures of our coarse-graining and multiscale modeling procedure for polymer melts.McCarty2 () By comparison with united atom simulations of polymer mixtures, our approach shows good agreement between the two descriptions. Furthermore, our procedure has the advantage of being analytically solved, so that thermodynamic parameters defining the free energy potential appear explicitly in the formalism, making the whole procedure general and transferible to other systems in different thermodynamic conditions. The key advantage of our method lies in the fact that it provides a quantitative picture of how the local details such as the specific chemical architecture and thermodynamic conditions are connected to the observed, global properties. While atomistic and united atom simulations are limited in the maximum length scale that can be simulated, our approach easily provides information on a wide range of length scales, once mesoscale simulations are combined with atomistic or monomer level simulations. The number of test systems presented in the paper is determined by the available united-atom simulations against which our procedure is tested; however, our procedure is directly applicable to other mixtures as well.

While the focus here is on the equilibrium structural properties of simple polyolefins, it is worth noting that many interesting systems could be investigated through our method and important problems approached. For example, the proposed multiscale procedure could be applied to study various aspects of binary blends, including critical behavior, mechanical properties and viscosity, the effects of non-random mixing, stiffness disparity, shape disparity, and different architectures, as well as the dynamics of demixing at equilibrium and under shear, interdiffusion in equilibrium and during phase separation, and the different propensity for crystallization in the dynamics of miscible blends approaching the phase transition. As one of the polymers in the mixture crystalizes it is important to investigate, for example, the kinetics of interface formation and structure. Additionally the procedure could be extended to encompass mixtures including block and homo-copolymers, blends close to surfaces and in thin films, and mixtures of colloids and polymers. Although the dynamics of coarse-grained systems is unrealistically accelerated, we have recently proposed a new first-principle method to rescale the dynamics of the coarse-grained system to its atomistic value, enabling dynamical properties of these systems to be investigated with quantitative precision.Lyubimov () This method can be useful in evaluating long-time dynamics directly from the mesoscale simulation.

The paper begins with a brief overview of the theory for coarse-graining of polymer blends, from which we calculate the effective potential, , input to the mesocale molecular dynamic simulations of the coarse-grained system. Mesoscale simulations of the coarse-grained polymer mixture, where chains are represented as soft colloidal particles, are presented in the following section. Those simulations provide the large-scale information, which is combined with local-scale united-atom simulations by means of our multiscale modeling procedure. Next, we present our results in the form of the total correlation function, with for each pair of components, which is predicted and compared with data of full united-atom simulations. First, we present correlation functions of the mixture in the limit where to demonstrate the tenor of our scale bridging approach. This is followed by an extension of the model to include finite temperature effects, where we present the change in the concentration structure factor upon inclusion of realistic values for the parameter. A brief discussion concludes the paper.

Ii Methodology

ii.1 Effective Pair Potentials for Polymer Mixtures

An analytical coarse graining procedure for binary polymer blends was derived by Yatsenko, et al.YAPRL (), to formally map the blend onto a mixture of soft-colloidal particles. Here, we briefly review the results of that work. The polymer species (A and B) in the blend are characterized by a radius of gyration, and , a total density of monomers, , and a segment length and , where and are the number of chain units. The volume fraction of the A component is given by and the ratio of segment lengths is given by . The main physical quantity of interest in the study of complex liquids, is the pair distribution function, or the related total correlation function, as from this quantity any thermodynamic property of interest can be directly calculated.McQuarrie () Here the total correlation function has contributions from three terms: , , , representing the three types on interaction of a mixture of particles A and B. Analytical expressions for the total correlation functions at the coarse grained level are obtained by solving a generalized Ornstein-Zernike (OZ) matrix relation


where H(k) is the matrix of total intermolecular partial correlation functions (pcfs), C(k) is the matrix of direct pcfs, and represents the matrix of intramolecular pcfs. Solving Equation 1 gives the relation for a generic pairYAJCP ()


where for each type of chain is the intramolecular form factor, and the superscripts, , , or , denotes center of mass-center of mass, monomer-center of mass, or monomer-mononomer interactions, respectively. The derivation of Equation 2 is an extension of the procedure outlined by Krakoviack, et al.K2002 () for homopolymers and formally connects the center of mass (cm) distribution functions to monomer-monomer (mm) distribution functions. Both the intramolecular form factors in Equation 2, namely the monomer-monomer and the monomer-cm , are expressed in analytical forms.Yamakawa (); DoiEd () Yatsenko, et al. introduced a new analytical form for the monomer total correlation functions, by extending the PRISM-blend thread model of Tang and SchweizerTANG () to include asymmetries in local chemical structure and flexibility.YAPRL ()

Following this procedure, analytical forms of , coarse-grained at the cm level, can be calculated from Equation 2. In real space they are given by:






where is the average correlation hole length scale, and is an average radius-of-gyration. Here, identifies the length scale for the density and concentration fluctuation correlations. The concentration fluctuation of a polymer blend is defined as


where we introduce the miscibility parameter for which the quantity is analogous to the Flory-Huggins interaction parameterYAJCP () such that . The concentration fluctuation, , diverges at the spinodal temperature where . The average segment length is , while the density fluctuation, , is defined as and . Although for compressible blends, there is, in principle, more than one parameter resulting from the distinct partial osmotic compressibilities of a two component blendSchweizer1995 (), in this work we adopt a single interaction parameter which determines the length scale of concentration fluctuations. The choice of a single adjustable parameter is appropriate in this case, since the aim is to obtain quantitative agreement with experiment, and the adoption of a single “effective” maintains a straightforward connection to SANS experiments. The notion of extracting a single effective parameter was also invoked Dudowicz et al. in a general lattice cluster theory analysis of compressible blends,Dudowicz () and by Schweizer in connecting the PRISM blend approach with methods used in experimental SANS analysis.SchweizerMacro () While we initially set , the extension to more realistic models is achieved by using the experimental parameter within this general framework, which has been shown to represent well both upper and lower critical phase diagrams.YAJCP () Results presented in Section III B of this paper are examples of both kinds of systems. In the athermal regime, however, the theory is presently being implemented to reproduce quantitatively and in a self-consistent way, the experimental values of workinprog ().

Finally, it should be noted that here we are only interested in capturing the global structure and thus the use of the simple analytic PRISM thread result is adequate as it accurately captures the structure of the polymer at distances greater the and correctly predicts the correlation hole. While the PRISM thread result does not capture local effects such as solvation shells, such detail is not relevant at the coarse-grained level since at this level, local interactions are averaged out. Local packing in the structural g(r) is introduced later using the multiscale modeling procedure discussed later, where local scale information from UA MD simulations is paired with global information obtained from course-grained simulations.

By adopting the hypernetted-chain (HNC) closure relation, the effective pair interaction potential, , is connected to the pcfs by the relationship


where is the direct pcf for a mixture of soft colloidal particles, defined in reciprocal space as


The chain density, , of chain type, A, is given as , and for chain type, B, as . For a binary mixture the static structure factors, , are given by

where the total chain density , and is the determinant of the mesoscopic static structure factor matrix. Whereas at the monomer level of description, molecular closures are required to ensure the correct scaling of the parameterSchweizerYethiraj (), the use of the site level HNC is appropriate here because at this level of coarse-graining, the polymers are treated as a simple liquid of soft-colloids.LouisPRL ()

By inserting Equations 3 and 9 into equation 8, the effective pair potentials are calculated, which are input to the mesocale simulations of the coarse grained binary mixture. The potential so obtained explicitly depends on the structural parameters of the polymer, such as , , , and . Since these parameters enter into the UA description, they do not represent additional parameters needed in moving to a higher level of coarse-graining. In other words, the potential is state dependent, being a free energy obtained from the monomer frame of reference; however, it is fully transferable to different systems since it is analytically derived in a well-defined manner. In addition to these structural parameters, there is the additional parameter, , that enters into the mesoscale description through Equation 7, and describes the monomer-specific interactions that drive phase separation. Once these parameters are defined, mesoscale simulations may be performed and structural properties calculated for any system of interest.

ii.2 Mesoscale Simulations of Coarse-Grained Polymer Mixtures

In continuing our multi-scaling approach, mesoscale simulations (MS) of various binary polymer mixtures with were performed using the effective potential calculated in the previous section. These simulations provide a reference system representative of the case in which the blend is well mixed and fluctuations are small. Systems investigated were blends of polyethylene (PE), polyisobutylene (PIB), and polypropylenes in their head-to-head (hhPP), isotactic (iPP), and syndiotactic (sPP) forms. Table LABEL:TB:1 shows the simulation parameters used in both UA and MS descriptions for each polyolefin blend studied.

Blend [A/B] [sites/] [] []
hhPP/sPP 0.50 0.0332 12.18 13.87 1.14
hhPP/PE 0.50 0.0332 12.32 16.48 1.34
PIB/PE 0.50 0.0343 9.76 16.38 1.68
PIB/sPP 0.50 0.0343 9.76 13.78 1.41
iPP/PE 0.75 0.0328 11.33 16.69 1.48
hhPP/PIB 0.50 0.0343 12.41 9.69 1.28
Table 1: Simulation Parameters for Polyolefin Blends (

Our MS MD simulations were performed on systems of point particles evolving in the microcanonical (, , ), ensemble. Initially all particles were placed on a cubic lattice with periodic boundary conditions, where the type of particle (/) occupying a particular lattice site was determined at random. The potential and its corresponding derivative were entered as tabulated inputs to the program, and linear interpolation was used to determine function values not supplied as the algorithm proceeds. Each site was given an initial velocity pooled from a Mersenne Twister random number generator,Matsumoto () and the system subsequently was evolved using a velocity Verlet integrator. We used reduced units such that all the units of length were scaled by () and energies were scaled by .

Equilibrium was induced by a quenching procedure, in which velocities were rescaled at regular intervals to maintain the desired temperature. Proper equilibration was verified by observing a Maxwell-Boltzmann distribution of velocities, a steady temperature, a stabilized Boltzmann H-theorem function, and a decayed translational order parameter. Once velocity rescaling was discontinued during the production stage, we monitored the fluctuating steady temperature to ensure the system remains in equilibrium throughout the duration of the simulation.

During the production stage, trajectories were collected over a traversal of . A typical MS MD simulation included particles, evolving for 30,000 computational steps. Several simulations were run using the LONI TeraGrid systemTeraGrid () to facilitate performing numerous simulations at a time. As a benchmark, a typical run on a single CPU workstation took h of wall clock time; however, since our codes are not yet subject to any parallelization process, we stress that this represents an underestimate of the potential gain in computational time as opposed to running full atomistic or UA simulations.

Blend [A/B] Particles [] []
hhPP/sPP 2744 199.07 166.61
hhPP/PE 5324 246.21 166.61
PIB/PE 4096 218.66 164.91
PIB/sPP 5488 245.19 164.91
iPP/PE 1728 168.57 167.27
hhPP/PIB 3456 230.73 164.91
Table 2: Mesoscale Simulation (MS-MD) Particle Number and Box Dimension Compared to UA Box Dimension. All UA simulations are for 1600 chains.

Mesoscale simulations provide the center-of-mass total correlation functions that describe the polymer mixtures on the large scale and are readily calculated from the simulation coordinates. As an exmple we show in Figure 1(a) the plot of for a 50:50 mixture of hhPP/sPP (). Data from mesoscale simulations and theoretical predictions are compared against united atom simulations for the center-of-mass total correlation functions and show an excellent agreement. Analytical theory, mesoscale simulations, and united atom simulations are all consistent in depicting the structure of the fluid on the length scale of the polymer radius of gyration and larger. Although all of the structural information is already contained in the analytic expression, the MS-MD simulations are useful as they can provide further information, for example, on the dynamics of the system close or far from equilibrium, both of which could be in principle investigated with our coarse-grained systems.

Figure 1: (a) Plot of for hhPP/sPP obtained from mesoscale simulation (open red circles). Comparison with UA data (filled circles) and theoretical predictions (solid line) shows quantitative agreement. (b) Plot of calculated using the inverse mapping procedure, Equation 11, (open red circles) compared to data from the full UA MD simulation (solid circle).

However, while describes the center of mass correlations between particles, when we are concerned with the structure of the liquid on the local scale, it is necessary to look at the monomer total correlation functions. Although the large scale information is completely provided by the mesoscale information, the information on the local-scale structure is averaged out during the coarse-graining. To take advantage of the gain in computational time given by mesoscale simulation, it is convenient to combine the large scale information from the latter with a short united atom simulation that describes the local structure. The two are merged together through a multiscale modeling procedure as presented in the next section. The advantage in performing this two-steps procedure is that the simulation with atomistic resolution can now be limited to short simulation runs and to a small ensemble of molecules, as the needed large-scale, long-time information comes from the fast mesoscale simulations.

ii.3 Recovering the Monomer Level Description: Implementation of the Multiscale Modeling Procedure for Polymer Blends

In our multiscale modeling we combine the mesoscale simulation, which correctly captures the large scale, global structure, with UA simulations, needed to determine the local properties of the polymer. By this procedure we aim at bridging the different length scales characteristic of the mixture and obtaining a complete description of the polymer blend structure at a defined temperature and composition.

The key to this procedure lies in extracting the monomer-monomer correlation function, from the coarse-grained correlation function, , by inverse mapping, using the reciprocal relation:


Here, is obtained from the mesoscale simulation, while the intramolecular form factors, , entering Equation 11 are calculated directly from short UA MD simulations. Since UA-MD simulations are needed only to obtain the intramolecular distributions, these simulations need only be for system box sizes on the order of the radius of gyration, which is a significant advantage over obtaining the global structure from UA-MD alone which requires much larger box sizes to accurately calculate the intermolecular distributions. The calculated is valid only for small k values, where the coarse-grained description applies, and begins to diverge as approaches zero at large k. This is shown in Figure 1(b), where we compare directly calculated using Equation 11 with UA-MD data of the same system.

To include local scale information, the total correlation function at small , obtained following the described procedure, is combined with from UA simulation for large values of . The key point in the procedure is to define a cut-off length, , at which the two descriptions are to be combined.

Because the local scale information is averaged out in the coarse-grained description, it is important to ensure that the whole intramolecular description is accounted for by the UA simulation. To estimate the length scale at which local intramolecular contributions become negligible, we evaluate the fraction of intra to total site/site contacts. For a given component of type in the mixture, this function is defined as


where is the number of type intramolecular contact sites


while the total number of site/site contacts is given by


In our procedure, the integral in Equation 13 is computed numerically using UA data for .

Selecting the correct cut-off distance is an important step of the procedure. Choosing a larger cutoff radius , smaller , includes more information from UA simulations, increasing the precision of the calculation, but it increases the computational time, partially defeating the purpose of adopting a multiscale procedure. In the case of polymer melts we observed that a value of gives good precision in the calculated when compared with full united atom simulation data, while retaining a reasonable speeding up of the calculations.McCarty2 () For polymer blends we adopt the same criterium and we evaluate a posteriori the agreement obtained for the blend samples analyzed in this paper. In Table LABEL:TB:NtBlnd we report, for both AA and BB type interactions, the effective radius, , for combining simulations, which lies in the intermediate length scale on the order of a few units and is unique for each blend. Also reported is the total number of site/site contacts within the radius determined by the value , evaluated using Equation 14 at . Dividing by the number of united atoms per chain (96 in our case) gives the number of chains, , needed in a UA simulation to produce the statistical information necessary to obtain the local structure. This is important since it determines how “short” the UA-MD simulations must be without losing pertinent information about local correlations.

Large-scale and local-scale information were combined at the chosen , and the procedure was repeated for each of the mixtures shown in Table LABEL:TB:1. The correlation coefficient between determined from our multiscaling procedure and from full UA simulations is calculated for the case when the multiscale simulations are combined at . These values are also presented in Table LABEL:TB:NtBlnd, showing that once is defined, the multiscale procedure provides an accurate way of obtaining the total correlation function for binary blends while circumventing the need for prohibitively large atomistic simulations. For the case of AB type interactions, an average radius between AA and BB data is used to determine where the simulations should be combined. While this approach allows one to explore a large range in the degree of polymerization, the method becomes impractical for liquids of long chains, for which intramolecular effects will remain long-ranged, since the cutoff length scale occurs on the order of . Moreover, for long chains, entanglement effects enter the dynamics and have to be accounted for. Therefore, for large N systems it is advisable to include a third, intermediate, level of coarse-graining, with a cut-off of the order of (and larger than) the persistence length and the entanglement length scale.

Iii Results

iii.1 Total Pair Correlation Functions of Athermal Reference Blends

The method just discussed gives a systematic way of merging simulations to optimize the tradeoff between the gain in accuracy due to inclusion of UA simulation data and the gain in efficiency due to the coarse grained mesoscopic picture. This procedure works well as it yields total correlation functions in excellent agreement with UA simulations at a reduced computational cost.

System Type Corr. Coeff
hhPP/sPP AA 29.1 3427 36 0.9999
hhPP/sPP BB 28.7 3288 34 0.9999
hhPP/PE96 AA 28.7 3288 34 0.9998
hhPP/PE96 BB 27.6 2924 31 0.9999
PIB/PE96 AA 28.8 3432 36 0.9999
PIB/PE96 BB 27.2 2891 30 0.9999
PIB/sPP AA 29.1 3536 37 0.9995
PIB/sPP BB 29.6 3721 39 0.9999
iPP/PE96 AA 29.3 3458 36 0.9999
iPP/PE96 BB 27.6 2891 30 0.9998
Table 3: Effective Radius for Combining Simulations and the Number of Total Site/Site Contacts Evaluated at

Results from the multiscale procedure are reported in Figures 2 - 6. In Figure 2, the left panels show the total correlation function, , obtained by combining UA simulations for the short length scales (large k) with mesoscale simulations for the long length scales (small k), for the first mixture in Table LABEL:TB:1, hhPP/sPP. The vertical dashed line in the left panels represents the cut-off length scale at which the two simulations are combined, which corresponds to . The insert in the left panels depicts the peak representing the local structure, which is accurately determined by UA simulation. This peak depends on the geometry of the monomeric structure and changes as a function of the type of polymers involved in the mixture.

The right panels in Figure 2 show , i.e. the Fourier transform of the corresponding total correlation function, . As in the case for melts, we adopted a sampling step of in reciprocal space, as this is of the same order as the discontinuity in , which results from joining the two simulations. Because of the amplitude of the chosen sampling step, there is no effect on the Fourier transformed function . The pair distribution function calculated from the multiscale modeling procedure is identical to the coarse-grained analytical expression in the large-scale regime, but also includes the local scale solvation shells, which come from UA-MD and are not captured by the PRISM thread expressions. Both pieces of information are needed to provide the complete structural and thermodynamic picture of the system.

Figure 2 illustrates the spirit of our multi-scale modeling approach in which independent simulations representing the same system at different levels of coarse-graining can be combined to provide a complete description of the polymer. Analogous plots for the hhPP/PE and the PIB/PE mixture are reported in Figure 3 and Figure 4, respectively. Finally, results for the two mixtures of PIB/sPP and iPP/PE are presented in Figures 5 and 6. All systems in real space show quantitative agreement with UA-MD data but are obtained at a much more efficient computational time than running the full UA MD simulation. Furthermore, the procedure to obtain the pair correlation functions is entirely analytical, hence we do not utilize any optimization procedure or numerical fitting scheme to obtain consistency between the two descriptions, and thus the method is fully transferable. Figures 2-6 demonstrate the versatility of the approach for studying mixtures of polymers with subtly different chemical architecture, since the same multiscale procedure is applied in all cases.

Figure 2: (Left) Multiscale modeling: The left panels show the total correlation function, , for AA (top), BB (middle), and AB (bottom) interactions, for a mixture of 50:50 hhPP/sPP. The data over the range of small k values was obtained by mesoscale simulation, whereas over the large k range it was obtained by UA MD simulation. The inset depicts the local structure. The dashed line indicates the value at which the two simulations were combined. (Right) The correlation function, , after Fourier transform (solid red line) is compared with results from the full UA MD simulation (open symbols).
Figure 3: Same as in Figure 2 for a blend of 50:50 hhPP/PE.
Figure 4: Same as in Figure 2 for a blend of 50:50 PIB/PE.
Figure 5: Same as in Figure 2 for a blend of 50:50 PIB/sPP.
Figure 6: Same as in Figure 2 for a blend of 75:25 iPP/PE.

iii.2 Applications to Thermal Blends with Realistic Parameters

In the multiscale modeling procedure presented above, we set in our soft-colloidal representation, and the UA simulations to which we combine our mesoscale simulations are assumed to be well-mixed. This provides an efficient means of obtaining pair correlation functions for polymer mixtures close to the athermal limit. In this section we examine the assumption that the blends are well-mixed and discuss the application of our multiscale modeling procedure to modeling thermal mixtures where subtle local features such as monomer shape, branching, and site energetic asymmetries gives rise to non-trivial divergent fluctuations.

Blends of polypropylenes of different tacticities have been studied by Woo et al., and the parameter of blends of iPP/sPP and aPP/sPP was found to be nearly zero over the temperature range T=423-453 K.Woo () Thus, for the hhPP/sPP blend of Figure 2 it is reasonable to assume the effective parameter will be small, which is supported by the overall quantitative agreement in Figure 1 and 3. For the other systems investigated in this publication, many of the pertinent parameters may be determined from the available literature on polyolefin blends which is presented in Table 4.

Blend [A/B] () hhPP/PE 11footnotemark: 1 0.077 PIB/PE 22footnotemark: 2 0.163 iPP/PE 0.0133footnotemark: 3 0.360 hhPP/PIB 44footnotemark: 4 0.018

Reference 26 Reference 27 Reference 25 Reference 23

Table 4: Temperature Dependence of Polyolefin Blends

To demonstrate the extension of our approach to thermal mixtures, we use the literature values for the parameter evaluated at the simulation temperature, and perform MS-MD simulations for the mixtures: hhPP/PE, PIB/PE, and hhPP/PIB. The hhPP/PIB blend is particularly interesting because it exhibits a lower critical solution temperature (LCST) phase diagram, which can in principle be modeled using our multiscale approach as well.McCarty1 () Since the parameters used in our mutiscale modeling procedure are determined using a UA description where each group is represented as a effective site, care must be taken when calculating the effective parameter from Table 4 where is defined on a monomer basis, for which several groups are grouped as one monomer unit. Since is proportional to the free energy of mixing per site, the parameter must be normalized by the average number of united atom sites per monomerJARAMILLO () as defined in the literature (6 for hhPP/PEJEON (), 4 for PIB/PELEE (), and 4.8 for hhPP/PIBJARAMILLO ()). For hhPP/PE this corresponds to a value of ; for PIB/PE a value of ; and hhPP/PIB a value of . For the temperature range at which the united atom simulations were performed, the magnitude of in all cases is small, supporting our initial modeling of in the previous section. In Table 4 we present the thermodynamically relevant parameter, at , which provides an indication of how far the mixture is from the spinodal.

We calculated the total correlation function for the thermal mixtures of hhPP/PE, PIB/PE, and hhPP/PIB, using the multiscale modeling approach described above. These results are presented in the supplemental materialSupplemental () and look nearly identical to the athermal results shown in Figures 3 and 4. Although the correlation functions are nearly identical with the athermal limit because the blends are well mixed, subtle differences in the structure as a result of increased concentration fluctuations should be manifest in the distributions. To assess the changes in structure that result from increased fluctuation, the monomer level partial static structure factors can be calculated from the total correlation function, ,


where the monomer form factors were determined from UA simulations as in Equation 11 above. The structure factor measuring correlations in the relative concentration, , which diverges as the mixture approaches the spinodal, is expressed as a linear combination of these partial structure factors,


In small angle neutron scattering (SANS) experiments, the parameter is determined from fitting the partial structure factor to the random phase approximation (RPA) equation of de GennesdeGennes ()


where for convenience the monomer site volumes were set equal to one.

Figure 7: (A) The concentration fluctuation structure factor obtained from the multiscale modeling procedure (open circles) for hhPP/PE () at . For comparison the RPA equation evaluated at (solid line) and (dashed line) is shown. Filled circles represent the structure factor from the full UA simulation. (B) The same as part (A), except for PIB/PE (), for which the RPA equation, evaluated at (solid line) and (dashed line), is shown. (C) the same except for the hhPP/PIB blend, and the RPA equation is evaluated at (solid line) and (dashed line).

Figure 7 presents for the three different thermal blends obtained from our multiscale modeling procedure. The static structure factor calculated in this manner exhibits good agreement when compared to the RPA equation, which was evaluated at for hhPP/PE, for PIB/PE, and for hhPP/PIB, using the intramolecular form factors from UA MD simulations. The relevance of this comparison to the RPA formula, which is known to fit low wavevector scattering curves well, is that there is clearly a better fit using the experimental parameter than with the case, which demonstrates the consistency of the description and that our approach is able to capture the fluctuations in concentration that arise in thermal polymer mixtures even at the relatively high temperatures of these simulations. Also, the advantage of a multiscale approach is exhibited by Figure 7 since the low wave vector regime is determined by mesoscale simulations so that UA simulations only need to be performed on systems at length scales up to the cut-off length (dashed line). This is important since only the initial stages of the divergent behavior need to be captured by united atom simulations, thus the need for prohibitively large simulation boxes is circumvented. For comparison, , calculated for the full UA MD simulation is also shown (offset for clarity) and agrees with our multiscale results, demonstrating the consistency between the two approaches.

When experimental values of the parameter are unavailable, a commonly used alternative is to estimate the interaction parameter from the solubility parameters calculated from a pure melt. These may be obtained from the cohesive energy determined from MD simulations of the individual melts. For the systems investigated here, we used united atom simulations of pure melts provided by Jaramillo et al.JARAMILLO () The cohesive energy density can be calculated from the radial distribution functions for each UA site type, , where designate a particular C, CH, CH, or CH group, according to


where is the potential between nonbonded sites, for which we employ a Lennard-Jones potential with the corresponding TraPPE-UA parameters of Martin and SiepmannSiepmann (). In Equation 18 the integration is carried out over the attractive branch of the potential. The solubility parameter (), calculated for each of the polyolefins in this study is presented in Table LABEL:TB:Sol. For comparison we also present results obtained by Pütz, et al.Putz () from UA simulations of chains with 24 backbone carbons.

Polymer (MPa)
PE 14.4 13.7 ()
hhPP 13.4 13.3
iPP 12.5 12.4
PIB 13.9 14.1
sPP 13.0 12.2
Table 5: Solubility Parameter for UA simulations of polyolefin melts (N=96) at . Values for are computed from Ref. Putz () for chains with 24 backbone carbons.

From these solubility parameters, the Flory-Huggins parameter may be estimated according to


where the densities, , are the pure melt densities of species A or B. Table LABEL:TB:chi shows a comparison between the Flory-Huggins parameter determined from this conventional solubility parameters approach and the observed experimental values. A major limitation of this solubility approach is that the value of is always positive and UCST phase behavior is always predicted. Furthermore, because non-combinatorial entropy is ignored, the concentration dependence of cannot be predicted. Despite these assumptions, the solubility approach has been demonstrated to hold up reasonably well for a large variety of polyolefin blendsLohse (). For comparison, in Table LABEL:TB:chi we also present the parameter computed directly from the full UA simulation of a polymer blend, by an extrapolation of the three lowest wave vector points to on an Ornstein-Zernike plot (1/S(k) vs. ), analogous to that obtained in typical SANS experiments. The lowest k vector is determined by the size of the simulation box and in this case is Å. Due to this limitation, the effective parameter can be shifted by about while still maintaining a reasonable fit to the blend structure factor.HEINE ()

Blend [A/B]
hhPP/PE 0.0016 0.0048 0.0022
PIB/PE 0.0034 0.0012 0.0027
hhPP/PIB 0.00038 0.0012 0.0017
Table 6: Comparison between experiment, solubility parameter, and RPA approaches to determine the parameter ( for the polyolefin blends in Figure 7

Finally we note that much work has been done using the PRISM theory of Schweizer and Curro to determine both the qualitative and quantitative miscibility behavior of polymer blends.Schweizer1994 (); Schweizer1993 (); Schweizer1997 () For polyolefin blends of the type investigated here, realistic models utilizing PRISM expressions for the pair distribution function along with UA potential parameters have been shown to give reasonable estimates of for various polypropylene and polyethylene blends.Rajasekaran (); Clancy (); HEINE (); however this approach is not pursued further here since the purpose of the current publication is to demonstrate the implementation of the multiscale modeling procedure once the parameter is obtained.

Iv Conclusion and Outlook

In this paper we present a multiscale modeling procedure to simulate mixtures of polymer chains at large scales. The method is completely general and applies to mixtures of polymers with different molecular structures, at different thermodynamic conditions of temperature, mixture composition, and chain length. The procedure combines large-scale information from a mesoscale simulation with a more detailed model which captures the local structure of the blend. For the detailed model, we use a united atom representation which captures the local polymer structure. Our mesoscale model represents the blend as a mixture of soft colloidal particles and accurately describes the large-scale structure as exhibited in the center of mass radial distribution function. Once the mesoscale simulations are performed, the relevant chemical details are reinserted by implementation of the Ornstein-Zernike formalism. We then combine this information with short united atom simulations to obtain a complete description over all length scales of interest. In this article we test the approach by applying our method to several different polyolefin mixtures. We check the validity of our procedure by direct comparison of the pair correlation function against full united atom molecular dynamic simulations. Because the UA-MD simulations available for comparison are limited to high temperatures where the samples can be more easily equilibrated, the test of our procedure includes only high temperature samples. However our mesoscale simulations are easily performed at temperatures approaching the demixing temperature. Furthermore, the multiscale modeling procedure provides a straightforward method to circumvent large simulations when modeling thermal mixtures, which has the potential to extend the range currently available to computer simulations of polymers. At low temperatures, where equilibration of even the local UA-MD simulation may be difficult, our two step procedure could be implemented “on the fly” where information is transferred from the local UA-MD (which provides the packing information) to the global scale MS-MD and vice versa until consistency between the two is reached.

One advantage of the proposed method is that no optimized parameterization scheme is required to maintain consistency between the different levels of description. The approach is thus transferable to different systems and different conditions, while remaining quantitative enough to distinguish even subtle differences in chemical structure. Considering the advantages of the proposed multiscale approach, it should provide an important tool for investigating the dynamics of large-scale phenomena via MD simulation. Future work should focus on performing simulations of large systems approaching the spinodal where full UA MD simulations are prohibitively difficult.

V Acknowledgements

We acknowledge support from the National Science Foundation through grant DMR-0804145. Computational resources were provided by LONI through the TeraGrid project supported by NSF.


  • (1) L. H. Sperling Polymeric Multicomponent Materials: An Introduction (John Wiley and Sons, Inc., New York, NY, 1997).
  • (2) D. Frenkel, and B. Smith Understanding Molecular Simulation. From Algorithms to Applications (Academic Press, London, 2002).
  • (3) K. Binder (ed.) Monte Carlo and Molecular Dynamics Simulations in Polymer Science (Oxford University Press, New York, 1995).
  • (4) E. Jaramillo, D. T. Wu, G. S. Grest, and J. G. Curro, J. Chem. Phys. 120, 8883 (2004).
  • (5) D. Heine, D. T. Wu, J. G. Curro, and G. S. Grest, J. Chem. Phys. 118, 914 (2003).
  • (6) M. P. Allen and D. J. Tildesley Computer Simulation of Liquids (Oxford Science Publications: Oxford, U.K. 1992.)
  • (7) F. Muller-Plathe Chem. Phys. Chem. 3, 754 (2002).
  • (8) J. I. Siepmann, S. Karaboni, and B. Smith, Nature 365, 330 (1993); P. Padilla and S. Toxvaerd, J. Chem. Phys. 95, 509 (1991); M. Mondello and G. S. Grest, J. Chem. Phys. 103, 7156 (1995); M.Mondello, G. S. Grest, A. R. Garcia, and B. G. Sibernagel, J. Chem. Phys. 105, 5208 (1996).
  • (9) M. Praprotnik, L. D. Site, and K. Kremer, Annu. Rev. Phys. Chem. 59, 545 (2008).
  • (10) S. A. Baeurle J. Math. Chem. 49, 363 (2009).
  • (11) J. McCarty, I. Y. Lyubimov, and M. G. Guenza, Macromol. 43, 3964 (2010).
  • (12) R. L. C. Akkermans, J. Chem. Phys. 128, 244904 (2008).
  • (13) T. D. Sewell, K. Ø. Rasmussen, D. Bedrov, G. D. Smith, and R. B. Thompson, J. Chem. Phys. 127, 144901 (2007).
  • (14) R. D. Groot and P. B. Warren, J. Chem. Phys. 107, 4423 (1997).
  • (15) R. D. Groot and K. L. Rabone, Biophys. J. 81, 725 (2001).
  • (16) K. Titievsky and G. C. Rutledge, J. Chem. Phys. 128, 124902 (2008).
  • (17) I. Y. Lyubimov, J. McCarty, A. Clark, and M. G. Guenza, J. Chem. Phys. 132, 224903 (2010).
  • (18) J. McCarty, I. Y. Lyubimov, and M. G. Guenza, J. Phys. Chem B 113, 11876 (2009).
  • (19) G. Yatsenko, E. J. Sambriski, M. A. Nemirovskaya, and M. Guenza, Phys. Rev. Lett. 93, 257803 (2004).
  • (20) D. A. McQuarrie Statistical Mechanics (University Science Books, Sunsalito, CA, 2000).
  • (21) G. Yatsenko, E. J. Sambriski, and M.G. Guenza, J. Chem. Phys. 122, 054907 (2005).
  • (22) V. Krakoviack, J.-P. Hansen, and A. A. Louis, Europhys. Lett. 58, 53 (2002).
  • (23) H. Yamakawa, Modern Theory of Polymer Solutions (Harper and Row, New York, 1971).
  • (24) M. Doi, and S. F. Edwards, The Theory of Polymer Dynamics (Oxford University, Oxford, 1986).
  • (25) H. Tang and K. S. Schweizer, J. Chem. Phys. 105, 799 (1996).
  • (26) C. Singh and K. S. Schweizer, J. Chem. Phys. 103, 5814 (1995).
  • (27) J. Dudowicz, M. S. Freed, and K. F. Freed, Macromol. 24, 5096 (1991).
  • (28) K. S. Schweizer, Macromol. 26, 6033 (1993).
  • (29) Manuscript in preparation.
  • (30) K. S. Schweizer and A. Yethiraj, J. Chem. Phys. 98, 9053 (1993); A. Yethiraj and K. S. Schweizer, J. Chem. Phys. 98, 9080 (1993).
  • (31) A. A. Louis, P. G. Bolhuis, J. P. Hansen, and E. J. Neijer, Phys. Rev. Lett. 85, 2522 (2000).
  • (32) M. Matsumoto and T. Nishimura, ACM Trans. on Modeling and Computer Simulation 8, 1, 3-30 (1998).
  • (33) C. Catlett, et al. HPC and Grids in Action; L. Grandinetti, Eds. IOS Press ‘Advances in Parallel Computing’ series; Amsterdam, 2007.
  • (34) E. M. Woo, K. Y. Cheng, Y-F. Chen, C. C. Su, Polymer 48 5753 (2007).
  • (35) H. S. Jeon, J. H. Lee, and N. P. Balsara, Macromol. 31, 3328 (1998).
  • (36) J. H. Lee, N. P. Balsara, A. K. Chakraborty, R. Krishnamoorti, and B. Hammouda, Macromol. 35, 7748 (2002).
  • (37) See Supplemental Material Document No. for additional figures.
  • (38) P. G. de Gennes, Scaling Concepts in Polymer Physics (Cornell University Press, Ithaca, NY, 1979).
  • (39) M. G. Martin and J. I. Siepmann, J. Phys. Chem. B, 103, 4508 (1999).
  • (40) M. Pütz, J. G. Curro, G. S. Grest, J. Chem. Phys. 114, 2847 (2001).
  • (41) R. Krishnamoorti, W. W. Graessley, G. T. Dee, D. J. Walsh, L. J. Fetters, and D. J. Lohse, Macromol. 29, 367 (1996).
  • (42) K. S. Schweizer and J. G. Curro, Adv. Polym. Sci. 116, 319 (1994).
  • (43) K. S. Schweizer, Macromol. 26, 6050 (1993).
  • (44) C. Singh and K. S. Schweizer, Macromol. 30, 1490 (1997).
  • (45) J. J. Rajasekaran, J. G. Curro, J. D. Honeycutt, Macromol. 28, 6843 (1995).
  • (46) T. C. Clancy, M. Pütz, J. D. Weinhold, J. G. Curro, and W. L. Mattice, Macromol. 33, 9452 (2000).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description