- Email: [email protected]
Multiscale modeling of polystyrene in various environments
Available online at www.sciencedirect.com
Fluid Phase Equilibria 261 (2007) 35–40
Multiscale modeling of polystyrene in various environments Qi Sun a , Florence R. Pon a,b , Roland Faller a,∗ a
Department of Chemical Engineering & Materials Science UC Davis, Davis, CA, United States b Intel Corporation, Folsom, CA, United States Received 31 March 2007; received in revised form 9 May 2007; accepted 11 May 2007 Available online 24 May 2007
Abstract We develop a systematic multiscale modeling approach for understanding the dynamics of polystyrene in the melt and blended with polyisoprene. The dynamics of polystyrene chains at chain length of 15 monomer in both environments are investigated using both atomistic and mesoscale models. The dynamics of polystyrene with chain length of 30 monomers and longer are illustrated by mesoscale models. The iterative Boltzmann inversion (IBI) method is used to derive the mesoscale models based on the atomistic simulations. The polystyrene chains clearly show dynamical heterogeneity along the chain and among chains in both models. The reorientation of PS chains speeds up with the increase of PI in the blends. We observe the onset of phase separation in atomistic simulations of polyisoprene–polystyrene blends at chain length of 15 monomer, and find full phase separation using mesoscale models for blends with chain length of 30 and above. We show that phase separation cannot be achieved with simple mixing rules. The phase separation morphologies could be cylindrical or spherical depending on the concentration ratio of the two components. The densities of the blends and their temperature dependencies are characterized. © 2007 Elsevier B.V. All rights reserved. Keywords: Polystyrene; Polyisoprene; Multiscale modeling
1. Introduction Atomistic simulations are required to understand the behaviors of polymers with precise details. However, there are many phenomena happening on larger scales where the required computer time in an atomistic model is prohibitively long. This applies, e.g., to entanglement effects [1] or to phase separation [2,3]. As chain lengths increase, polymer dynamics can only be captured by a mesoscale model. There are extensive efforts in mapping polymer systems from atomistic simulations onto the mesoscale [4–8]. The iterative Boltzmann inversion (IBI) [9] is one of the mapping techniques to derive mainly non-bonded potentials in a systematic way. It leads to models more suitable than simple bead-spring models [10], since it still keeps the chemical identity of polymer chains while speeding up the simulations. In this paper we discuss the comparison of atomistic and mesoscale models for polystyrene in the melt and the blend. The entanglement effects of polymers in simulations have been
∗
Corresponding author. E-mail address: [email protected] (R. Faller).
0378-3812/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2007.05.020
investigated in detail [11–13]. We first show some data of PS homo-polymers before we come to the main focus of the paper which is PS in the blend with polyisoprene (PI). PI–PS blends have been known to phase separate experimentally [14,15]. Our research provides a systematic approach in answering the above questions. 2. Simulation methodology 2.1. Atomistic simulation 2.1.1. Force field All atomistic simulations are conducted at 450 K and ambient pressure with the GROMACS package; we used a 20 ns equilibration time. The polyisoprene (PI) model has been taken from previous simulations [16]. Only, here the bond lengths are constrained using the LINCS algorithm [17] and we simulate cis-PI. For the force field of polystyrene (PS), we use the model described in [18]. All bonds are constrained, and all 1–4 non-bonded interaction are excluded. Constant temperature and pressure are enforced using the weak coupling method [19] with time constants of 0.2 ps for temperature and 1.0 ps for pressure. The time step is 2 fs and configurations are saved every 2000
36
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
Table 1 Polymer systems of mesoscale simulations PS melt
PI–PS blends
Chain length (-mer)
Number of chains
Systems
Weight ratio
Number of chains
7 10 15 30 40 50 70 80 90 100 120 160 240
48 48 48 48 48 48 48 48 48 48 24 24 12
7PI7–48PS7 72PI10–48PS10 72PI14–48PS15 72PI15–48PS15 36PI30–24PS30 36PI30–24PS45 36PI30–24PS60 36PI30–24PS80 36PI45–24PS30 36PI60–24PS30 36PI80–24PS30 36PI60–24PS60 36PI80–24PS80
50:50 50:50 50:50 50:50 50:50 40:60 34:66 28:72 60:40 66:34 72:28 50:50 50:50
120 120 120 120 60 60 60 60 60 60 60 60 60
steps. All chains are of length 15 monomers, PS chains are atactic with a 50% probability of racemo and meso dyads. The number of chains for each simulation is listed in Table 1. Density, end-to-end distance and radius of gyration are used to test for equilibration. 2.2. Mapping Mapping is the procedure to develop a mesoscale model from an atomistic parent simulation. Throughout the mapping process, mesoscale models are set up with a topology of simple bead-spring models, i.e., we have one interaction site per monomer and simple connections by bonds. All bonds, angles, and associated forces of the mesoscale models are derived from the atomistic model. Details can be found in references [1,2,20,21]. The interaction centers of the PS monomers are located at the backbone carbons connecting to the phenyl groups. The interaction centers of PI are located at the centers of the single bonds connecting two neighboring atomistic monomers [21]. Deriving non-bonded effective potentials is the most demanding part in contrast to the bond and angle parameters. The non-bonded potentials are optimized against the atomistic simulation until the radial distribution function generated from the mesoscale model is consistent with the atomistic simulation. The Iterative Boltzmann Inversion technique [9] is used to generate meaningful potential parameters. The mapping for polymer blends is more demanding than for pure polymer melts in that we need to reproduce the distribution functions of three pairs: PI–PI, PS–PS, PI–PS at the same time compared to only the PS–PS pair. Details of the iterations can be found in references [2,20]. 2.3. Mesoscale simulation PS melts and PI–PS blends are investigated using the MDSPHERICAL package [22] with a fixed density of 5.7499 and 6.4175 monomers/σ 3 . The equation of motion is integrated with a time step of 0.005τ. The bond lengths of PI and PS are 0.2345
and 0.2545 nm, respectively. The cutoff for the non-bonded interaction is 1.910σ (which is the same as has been used during the optimization). The non-bonded numerical potentials of the 36PI–24PS at chain length of 15 monomer have been applied to all the other blend systems. We additionally performed some simulations where we used empirical mixing rules for the PI–PS interaction. We used two different scenarios. First we applied an algebraic average of the PI–PI and PS–PS interaction for the unlike PI–PS interactions. Furthermore we tested a geometrical average in spirit of Lorentz–Berthelot mixing rules [23]. At large distances where the PS–PS interaction is slightly negative we set the non-bonded potentials to zero. These simulations were set up with equilibrated structures from the systematically mapped simulations. Various systems of PS melts and blends have been simulated and are listed in Table 1. We use atactic PS and 1,4-cisPI. 2PS–58PI, 10PS–46PI, 36PI–24PS, 37PS–19PI, 2PI–36PS, 48PS which corresponds to a PS concentration ranging from 5, 25, 50, 75, 95 up to 100%. In our nomenclature, e.g., 2PS–58PI stands for a system with 2PS and 58PI chains; 36PI30–24PS60 for a blend with 36PI chains at chain length of 30 monomer and 24PS chains at chain length of 60 monomer. 3. Results and discussion 3.1. Mapping Mapping is conducted based on atomistic simulations. It is not feasible to conduct corresponding atomistic scale simulations for longer chains. Only atomistic simulations of a pure PS melt with 48 chains at chain length 15 monomers and a 36PI15–24PI15 blend with a weight concentration ratio of 50:50 are set up to derive the non-bonded potentials for all mesoscale melts and blends. Through this mapping process, we keep the chemical structures while coarsening the atomistic scale model. Starting from the atomistic model of pure polystyrene oligomers, we investigate the local structure, particularly the radial distribution function, and the bond, angle distribution. The final RDF
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
Fig. 1. Diffusion constants at various chain lengths.
of the optimized mesoscale model is well consistent with the atomistic target. The discrepancy is within 5%. The generalization to binary blends requires three RDFs, that is, PI–PI, PS–PS and PI–PS pairs. RDF pairs are intermolecular, which means that pairs from the same polymer chain are excluded. Details of the iteration process and results can be found in [2,3]. The iteration process of the binary blend is more demanding compared to the melt as the different RDFs are interrelated. The maximum deviation for the non-bonded potentials throughout the three pairs is less than 10%. Note that for most of the RDF the deviation is well below this upper limit. 3.2. Dynamics of PS chains in the melt and blends The dynamics of polystyrene melts can now be investigated at various chain lengths ranging from 15 to 240 monomers. The crossover to entanglement behavior could be captured by several manifestations: mean-squared displacements, Rouse mode analysis and diffusion constants. A detailed discussion of these characteristics can be found in reference [1]. Here, we show the diffusion constant of the chain center of the mass as an example. With the onset of entanglements, the diffusion constants drops
37
from D ∝ N−1 to D ∝ N−2 [24]. Diffusion constants of various systems are calculated to verify the crossover from the unentangled to entangled dynamics (cf. Fig. 1). The entanglement length of this mesoscale model determined as around 85 monomers by the diffusion constants investigation. An experimental value of Ne ≈ 130 is determined by measuring the plateau modulus under oscillatory shear at a temperature of 413 K [25]. As our value of Ne ≈ 85 corresponds to atomistic parent simulations at 450 K and the entanglement effects decrease with the increase of temperature, we expect the discrepancy is minor. To understand the PS chain dynamics in blends, the reorientation correlation function of end-to-end PS vectors from both atomistic scale and mesoscale simulations are plotted in Fig. 2. It is very expensive to simulate the overall reorientation of polymer chains using an atomistic scale model atomistic scale. Therefore, the polymer system with chain length 15 monomers is selected for the comparison as it provides a reasonable balance between efficiency and minimizing end effects. There is clearly a strong dynamical heterogeneity among the total 24PS chains as a wide distribution of reorientation behaviors is observed. This has been seen in atomistic simulations earlier [26]. The reorientation functions are shown as a function of time using picoseconds (ps) in the atomistic scale and “unit time” τ in mesoscale. We convert the τ units to real time unit by calibrating the diffusion constants in both models. We have used the 36PI–24PS chain length of 15 since both scale models are available. The diffusion constant D of PS for this specific system in the atomistic scale is D = 0.018 × 10−5 cm2 /s, while in the mesoscale we obtain D = 3.378 × 10−16 cm2 /τ leading to τ = 2 ns is derived. This conversion confirms that mesoscale models can capture phenomena on much longer scales and speed up the simulations dramatically. We can classify the dynamics of PS chains according to reorientation correlation function based on the first order Legendre (P1 (t)) which is defined as (t) · u (0). P1 (t) = u
(1)
Reorientation of end-to-end vectors could be one of the ways to characterize the PS whole chain movement and it is used here.
Fig. 2. Orientation correlation functions of PS end-to-end vectors from atomistic scale (left) and mesoscale (right) simulations in the 36PI and 24PS systems. At long times we can expect significant statistical error.
38
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
superatom chain into two parts, the middle and end part. The middle group includes the central three atoms or superatoms based on atomistic or mesoscale simulations, while the end group is composed by the three monomers closest to the ends on either side. The corresponding correlation functions are averaged to understand the heterogeneity along the chain in the dynamics of segments. In both atomistic and mesoscale models, the ends of the chains relax much faster than the middle ones as expected (Fig. 4). 3.3. Phase separation in PI–PS blends
Fig. 3. Rotation correlation function of PI (left) and PS (right) end-to-end vectors at various polymer systems.
Computer simulations can address not only the average properties of the system, but also the distribution over any observable of interest. The PS chains could, e.g., be classified as “fast”, “intermediate” and “slow” by splitting the total 24 chains into groups according to their rotation speed as seen from Fig. 2. We now perform a dynamic analysis for the atomistic simulations of PI–PS blends. We select systems of 58PI–2PS, 48PS, 36PI–24PS as examples. By comparing the average reorientations of PI and PS from 36PI–24PS, we confirmed the experimental suggestion [27] that the correlation times of PS are much slower than those of PI as shown in Fig. 3. The dynamics of polymer chains are not only determined by their own identity, but also strongly affected by their surrounding chains. By investigating the relaxation process of PI in 58PI–2PS and 36PI–24PS or PS in 48PS and 36PI–24PS, we validate that polymer chains movements do speed up or slow down depending on the average dynamics of neighboring chains and as shown above this leads to a wide variety of individual re-orientation functions. Fig. 3 clearly shows that the PI chains loose their initial orientation much faster with the increase of the PI concentration and PS chains become slower with the increase of PS concentration. The segmental motion of PS chains has also been characterized. We divide the atomistic 15 monomer or mesoscale 15
Using our mesoscale models we can now embark on a study of the phase behavior. We observe clear lamellae pattern for an equiconcentrated (by mass) mixture for 36PI–24PS30. The systems take on a cylindrical shape as the concentration of monomers changes to 60:40 and 66:34 in the cases of 36PI45–24PS30 and 36PI60–24PS30, respectively. The morphology changes into a spherical shape as the concentration goes to 72:28 as in the system of 36PI80–24PS30 [2,3]. To characterize the lamellar morphologies we calculate the density profiles of systems with a weight concentration of 50:50 in binary blends from both atomistic and mesoscale simulations. The chain lengths are 15, 30, and 60 monomers. We observe clearly phase separation behaviors in all cases, the only difference is the thickness of the interface region. As indicated in our previous work, a 15 monomer chain is too short to be classified as a polymer, it is only oligomers [1]. Polymers, at chain length of 30 or above, have to be used to investigate the characteristics of polymer blends. Fig. 5(b and c) are manifestations of a lamellar morphology that has formed during the phase separation as in two dimensions the densities are flat and the concentrations varies in one dimension. 3.4. Empirical mixing rules Let us now discuss the inapplicability of empirical mixing rules for this system. As mentioned above, we performed test simulations using algebraically and geometrically aver-
Fig. 4. Rotation correlation function of PS segmental vectors in 36PI15mer-24PS15mer from atomistic scale (left) and mesoscale (right).
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
39
Fig. 5. Density profile of 36PI–24PS at (a) chain length of 15 monomer from atomistic scale, (b) chain length of 30 monomer from mesoscale and (c) chain length of 60 monomer from mesoscale.
aged homo-interactions (PI–PI, PS–PS) as approximations for the hetero-interaction (PIPS), i.e., we approximated VPI–PS = 0.5(VPI–PI + VPSPS ) or VPI–PS = (VPI–PI + VPSPS )0.5 . At large distances where the VPS–PS is slightly negative we set VPI–PS = 0. Fig. 6 shows final results of systems which were set
up in morphologies equilibrated by the systematically mapped potentials. The empirical mixing rules clearly underestimate the phase separation tendency and lead to mixed systems in contrast to experimental results and the systematically mapped simulations. These comparison results clearly validate that an
Fig. 6. Snapshots of mixed systems using empirical mixing rules after being started in a morphology equilibrated by the systematic potential (a) 36PI6024PS60 using a geometrical averaging started from lamellae and (b) 36PI4524PS30 using an algebraic averaging started from a cylinder.
40
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
independent optimization of the hetero-interactions is inevitable in capturing the phase separation behaviors. 4. Conclusion Polymer melts and blends of chain lengths ranging from 7 up to 240 monomers with blends of various concentrations have been investigated from both atomistic and mesoscale models. We applied IBI in mapping the mesoscale models from atomistic simulations. With these mesoscale models, we reproduce the entanglement behaviors for longer PS chains and quantify the entanglement length of PS melt. Overall, the dynamics of PS chains is slower than that of PI and it speed up with the number of surrounding PS chains. There are clearly dynamic heterogeneity along the chain and among the chains. Phase separation of PI–PS blends is observed at chain length of 30 monomer and up. The morphology takes on a lamellar shape at concentration of 50:50, cylindrical at 60:40 and 66:34, or leads to a spherical shape at 72:28. The density profiles of blends are also characterized to provide quantitative results. It is noteworthy that in one of our earlier publications [2] we preliminarily suggested, based solely on RDFs, a cylindrical morphology for chains of length 80. After a more detailed analysis it turns out that that system is lamellar as well. It is clear that in principle the mesoscale potentials depend on concentration. Therefore, as we use the potential developed at a certain concentration we would expect that the boundaries between the different morphologies are only tentative. But it is clear that this approach is superior to the application of mixing rules as they lead to qualitatively wrong mixing behavior. Acknowledgments Financial support from the US Department of Energy, Office of Advanced Scientific Computing, under grant DE-FG0203ER25568 is gratefully acknowledged. References [1] Q. Sun, R. Faller, Crossover from unentangled to entangled dynamics in a systematically coarse-grained polystyrene melt, Macromolecules 39 (2) (2006) 812–820. [2] Q. Sun, R. Faller, Systematic coarse-graining of a polymer blend: polyisoprene and polystyrene, J. Chem. Theor. Comput. 2 (3) (2006) 607–615. [3] Q. Sun, R. Faller, Phase separation in polyisoprene/polystyrene blends by a systematically coarse-grained model, J. Chem. Phys. 126 (14) (2007) 144908. [4] W. Tsch¨op, K. Kremer, J. Batoulis, T. B¨urger, O. Hahn, Simulation of polymer melts. I. Coarse-graining for polycarbonates, Acta Polym. 49 (1998) 61–74. [5] J. Baschnagel, K. Binder, P. Doruker, A.A. Gusev, O. Hahn, K. Kremer, W.L. Mattice, F. M¨uller-Plathe, M. Murat, W. Paul, S. Santos, U.W. Suter, V. Tries, Bridging the gap between atomistic and coarse-grained models of polymers: Status and perspectives, Adv. Polym. Sci. 152 (2000) 141–156.
[6] H. Meyer, O. Biermann, R. Faller, D. Reith, F. M¨uller-Plathe, Coarse graining of nonbonded interparticle potentials using automatic simplex optimization to fit structural properties, J. Chem. Phys. 113 (15) (2000) 6264–6275. [7] R. Faller, F. M¨uller-Plathe, Multi-scale modeling of poly(isoprene) melts, Polymer 43 (2) (2002) 621–628. [8] R. Faller, Coarse-grain modeling of polymers, in: K. Lipkowitz, T. Cundari (Eds.), Reviews in Computational Chemistry, vol. 23, Wiley-VCH, 2007, pp. 233–262 (Chapter 4). [9] D. Reith, M. P¨utz, F. M¨uller-Plathe, Deriving effective mesoscale potentials from atomistic simulations, J. Comput. Chem. 24 (2003) 1624–1636. [10] G.S. Grest, M.-D. Lacasse, K. Kremer, A.M. Gupta, Efficient continuum model for simulating polymer blends and copolymers, J. Chem. Phys. 105 (23) (1996) 10583–10594. [11] K. Kremer, G.S. Grest, Dynamics of entangled linear polymer melts: a molecular-dynamics simulation, J. Chem. Phys. 92 (8) (1990) 5057–5086. [12] V.A. Harmandaris, V.G. Mavrantzas, D.N. Theodorou, M. Kr¨oger, J. ¨ Ramirez, H.C. Ottinger, D. Vlassopoulos, Crossover from the rouse to the entangled polymer melt regime: signals from long, detailed atomisitic molecular dynamics simulations, supported by rheological experiments, Macromolecules 36 (2003) 1376–1387. [13] S. Le´on, N. van der Vegt, L.D. Site, K. Kremer, Bisphenol a polycarbonate: entanglement analysis from coarse-grained md simulations, Macromolecules 38 (2005) 8078–8092. [14] Y. He, T.R. Lutz, M.D. Ediger, M. Pitsikalis, N. Hadjichristidis, E.D. von Meerwall, Miscible polyisoprene/polystyrene blends: distinct segmental dynamics but homogeneous terminal dynamics, Macromolecules 38 (2005) 6216–6226. [15] R. Koningsveld, W.J. MacKnight, Liquid–liquid phase separation in multicomponent polymer systems XXVII. Determination of the pair interaction function for polymer blends, Polym. Int. 44 (3) (1997) 356–364. [16] R. Faller, F. Muller-Plathe, M. Doxastakis, D. Theodorou, Local structure and dynamics of trans-polyisoprene oligomers, Macromolecules 34 (5) (2001) 1436–1448. [17] B. Hess, H. Bekker, H.J.C. Berendsen, J.G.E.M. Fraaije, LINCS: a linear constraint solver for molecular simulations, J. Comput. Chem. 18 (12) (1997) 1463–1472. [18] Q. Sun, R. Faller, Molecular dynamics of a polymer in mixed solvent: atactic polystyrene in a mixture of cyclohexane and n,n-dimethylformamide, J. Phys. Chem. B 109 (33) (2005) 15714–15723. [19] H.J.C. Berendsen, J.P.M. Postma, V. Gunsteren, D.W.F.A.J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (8) (1984) 3684–3690. [20] Q. Sun, R. Faller, Systematic coarse-graining of atomistic models for simulation of polymeric systems, Comput. Chem. Eng. 29 (11–12) (2005) 2380–2385. [21] R. Faller, D. Reith, Properties of polyisoprene model building in the melt and in solution, Macromolecules 36 (14) (2003) 5406–5414. [22] S. Hendrik Meyer, Private correspondence, Unpublished data, 2005. [23] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press Inc., 1987. [24] G.R. Strobl, The Physics of Polymers: Concepts for Understanding their Structures and Behavior, second ed., Springer, 1997. [25] L.J. Fetters, D.J. Lohse, D. Richter, T.A. Witten, A. Zirkel, Reviews connection between polymer molecular weight, density, chain dimensions, and melt viscoelastic properties, Macromolecules 27 (17) (1994) 4639–4647. [26] R. Faller, Correlation of static and dynamic inhomogeneities in polymer mixtures: a computer simulation of polyisoprene and polystyrene, Macromolecules 37 (3) (2004) 1095–1101. [27] Y. He, T.R. Lutz, M.D. Ediger, Segmental and terminal dynamics in miscible polymer blends, J. Chem. Phys. 119 (18) (2003) 9956–9965.
Fluid Phase Equilibria 261 (2007) 35–40
Multiscale modeling of polystyrene in various environments Qi Sun a , Florence R. Pon a,b , Roland Faller a,∗ a
Department of Chemical Engineering & Materials Science UC Davis, Davis, CA, United States b Intel Corporation, Folsom, CA, United States Received 31 March 2007; received in revised form 9 May 2007; accepted 11 May 2007 Available online 24 May 2007
Abstract We develop a systematic multiscale modeling approach for understanding the dynamics of polystyrene in the melt and blended with polyisoprene. The dynamics of polystyrene chains at chain length of 15 monomer in both environments are investigated using both atomistic and mesoscale models. The dynamics of polystyrene with chain length of 30 monomers and longer are illustrated by mesoscale models. The iterative Boltzmann inversion (IBI) method is used to derive the mesoscale models based on the atomistic simulations. The polystyrene chains clearly show dynamical heterogeneity along the chain and among chains in both models. The reorientation of PS chains speeds up with the increase of PI in the blends. We observe the onset of phase separation in atomistic simulations of polyisoprene–polystyrene blends at chain length of 15 monomer, and find full phase separation using mesoscale models for blends with chain length of 30 and above. We show that phase separation cannot be achieved with simple mixing rules. The phase separation morphologies could be cylindrical or spherical depending on the concentration ratio of the two components. The densities of the blends and their temperature dependencies are characterized. © 2007 Elsevier B.V. All rights reserved. Keywords: Polystyrene; Polyisoprene; Multiscale modeling
1. Introduction Atomistic simulations are required to understand the behaviors of polymers with precise details. However, there are many phenomena happening on larger scales where the required computer time in an atomistic model is prohibitively long. This applies, e.g., to entanglement effects [1] or to phase separation [2,3]. As chain lengths increase, polymer dynamics can only be captured by a mesoscale model. There are extensive efforts in mapping polymer systems from atomistic simulations onto the mesoscale [4–8]. The iterative Boltzmann inversion (IBI) [9] is one of the mapping techniques to derive mainly non-bonded potentials in a systematic way. It leads to models more suitable than simple bead-spring models [10], since it still keeps the chemical identity of polymer chains while speeding up the simulations. In this paper we discuss the comparison of atomistic and mesoscale models for polystyrene in the melt and the blend. The entanglement effects of polymers in simulations have been
∗
Corresponding author. E-mail address: [email protected] (R. Faller).
0378-3812/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.fluid.2007.05.020
investigated in detail [11–13]. We first show some data of PS homo-polymers before we come to the main focus of the paper which is PS in the blend with polyisoprene (PI). PI–PS blends have been known to phase separate experimentally [14,15]. Our research provides a systematic approach in answering the above questions. 2. Simulation methodology 2.1. Atomistic simulation 2.1.1. Force field All atomistic simulations are conducted at 450 K and ambient pressure with the GROMACS package; we used a 20 ns equilibration time. The polyisoprene (PI) model has been taken from previous simulations [16]. Only, here the bond lengths are constrained using the LINCS algorithm [17] and we simulate cis-PI. For the force field of polystyrene (PS), we use the model described in [18]. All bonds are constrained, and all 1–4 non-bonded interaction are excluded. Constant temperature and pressure are enforced using the weak coupling method [19] with time constants of 0.2 ps for temperature and 1.0 ps for pressure. The time step is 2 fs and configurations are saved every 2000
36
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
Table 1 Polymer systems of mesoscale simulations PS melt
PI–PS blends
Chain length (-mer)
Number of chains
Systems
Weight ratio
Number of chains
7 10 15 30 40 50 70 80 90 100 120 160 240
48 48 48 48 48 48 48 48 48 48 24 24 12
7PI7–48PS7 72PI10–48PS10 72PI14–48PS15 72PI15–48PS15 36PI30–24PS30 36PI30–24PS45 36PI30–24PS60 36PI30–24PS80 36PI45–24PS30 36PI60–24PS30 36PI80–24PS30 36PI60–24PS60 36PI80–24PS80
50:50 50:50 50:50 50:50 50:50 40:60 34:66 28:72 60:40 66:34 72:28 50:50 50:50
120 120 120 120 60 60 60 60 60 60 60 60 60
steps. All chains are of length 15 monomers, PS chains are atactic with a 50% probability of racemo and meso dyads. The number of chains for each simulation is listed in Table 1. Density, end-to-end distance and radius of gyration are used to test for equilibration. 2.2. Mapping Mapping is the procedure to develop a mesoscale model from an atomistic parent simulation. Throughout the mapping process, mesoscale models are set up with a topology of simple bead-spring models, i.e., we have one interaction site per monomer and simple connections by bonds. All bonds, angles, and associated forces of the mesoscale models are derived from the atomistic model. Details can be found in references [1,2,20,21]. The interaction centers of the PS monomers are located at the backbone carbons connecting to the phenyl groups. The interaction centers of PI are located at the centers of the single bonds connecting two neighboring atomistic monomers [21]. Deriving non-bonded effective potentials is the most demanding part in contrast to the bond and angle parameters. The non-bonded potentials are optimized against the atomistic simulation until the radial distribution function generated from the mesoscale model is consistent with the atomistic simulation. The Iterative Boltzmann Inversion technique [9] is used to generate meaningful potential parameters. The mapping for polymer blends is more demanding than for pure polymer melts in that we need to reproduce the distribution functions of three pairs: PI–PI, PS–PS, PI–PS at the same time compared to only the PS–PS pair. Details of the iterations can be found in references [2,20]. 2.3. Mesoscale simulation PS melts and PI–PS blends are investigated using the MDSPHERICAL package [22] with a fixed density of 5.7499 and 6.4175 monomers/σ 3 . The equation of motion is integrated with a time step of 0.005τ. The bond lengths of PI and PS are 0.2345
and 0.2545 nm, respectively. The cutoff for the non-bonded interaction is 1.910σ (which is the same as has been used during the optimization). The non-bonded numerical potentials of the 36PI–24PS at chain length of 15 monomer have been applied to all the other blend systems. We additionally performed some simulations where we used empirical mixing rules for the PI–PS interaction. We used two different scenarios. First we applied an algebraic average of the PI–PI and PS–PS interaction for the unlike PI–PS interactions. Furthermore we tested a geometrical average in spirit of Lorentz–Berthelot mixing rules [23]. At large distances where the PS–PS interaction is slightly negative we set the non-bonded potentials to zero. These simulations were set up with equilibrated structures from the systematically mapped simulations. Various systems of PS melts and blends have been simulated and are listed in Table 1. We use atactic PS and 1,4-cisPI. 2PS–58PI, 10PS–46PI, 36PI–24PS, 37PS–19PI, 2PI–36PS, 48PS which corresponds to a PS concentration ranging from 5, 25, 50, 75, 95 up to 100%. In our nomenclature, e.g., 2PS–58PI stands for a system with 2PS and 58PI chains; 36PI30–24PS60 for a blend with 36PI chains at chain length of 30 monomer and 24PS chains at chain length of 60 monomer. 3. Results and discussion 3.1. Mapping Mapping is conducted based on atomistic simulations. It is not feasible to conduct corresponding atomistic scale simulations for longer chains. Only atomistic simulations of a pure PS melt with 48 chains at chain length 15 monomers and a 36PI15–24PI15 blend with a weight concentration ratio of 50:50 are set up to derive the non-bonded potentials for all mesoscale melts and blends. Through this mapping process, we keep the chemical structures while coarsening the atomistic scale model. Starting from the atomistic model of pure polystyrene oligomers, we investigate the local structure, particularly the radial distribution function, and the bond, angle distribution. The final RDF
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
Fig. 1. Diffusion constants at various chain lengths.
of the optimized mesoscale model is well consistent with the atomistic target. The discrepancy is within 5%. The generalization to binary blends requires three RDFs, that is, PI–PI, PS–PS and PI–PS pairs. RDF pairs are intermolecular, which means that pairs from the same polymer chain are excluded. Details of the iteration process and results can be found in [2,3]. The iteration process of the binary blend is more demanding compared to the melt as the different RDFs are interrelated. The maximum deviation for the non-bonded potentials throughout the three pairs is less than 10%. Note that for most of the RDF the deviation is well below this upper limit. 3.2. Dynamics of PS chains in the melt and blends The dynamics of polystyrene melts can now be investigated at various chain lengths ranging from 15 to 240 monomers. The crossover to entanglement behavior could be captured by several manifestations: mean-squared displacements, Rouse mode analysis and diffusion constants. A detailed discussion of these characteristics can be found in reference [1]. Here, we show the diffusion constant of the chain center of the mass as an example. With the onset of entanglements, the diffusion constants drops
37
from D ∝ N−1 to D ∝ N−2 [24]. Diffusion constants of various systems are calculated to verify the crossover from the unentangled to entangled dynamics (cf. Fig. 1). The entanglement length of this mesoscale model determined as around 85 monomers by the diffusion constants investigation. An experimental value of Ne ≈ 130 is determined by measuring the plateau modulus under oscillatory shear at a temperature of 413 K [25]. As our value of Ne ≈ 85 corresponds to atomistic parent simulations at 450 K and the entanglement effects decrease with the increase of temperature, we expect the discrepancy is minor. To understand the PS chain dynamics in blends, the reorientation correlation function of end-to-end PS vectors from both atomistic scale and mesoscale simulations are plotted in Fig. 2. It is very expensive to simulate the overall reorientation of polymer chains using an atomistic scale model atomistic scale. Therefore, the polymer system with chain length 15 monomers is selected for the comparison as it provides a reasonable balance between efficiency and minimizing end effects. There is clearly a strong dynamical heterogeneity among the total 24PS chains as a wide distribution of reorientation behaviors is observed. This has been seen in atomistic simulations earlier [26]. The reorientation functions are shown as a function of time using picoseconds (ps) in the atomistic scale and “unit time” τ in mesoscale. We convert the τ units to real time unit by calibrating the diffusion constants in both models. We have used the 36PI–24PS chain length of 15 since both scale models are available. The diffusion constant D of PS for this specific system in the atomistic scale is D = 0.018 × 10−5 cm2 /s, while in the mesoscale we obtain D = 3.378 × 10−16 cm2 /τ leading to τ = 2 ns is derived. This conversion confirms that mesoscale models can capture phenomena on much longer scales and speed up the simulations dramatically. We can classify the dynamics of PS chains according to reorientation correlation function based on the first order Legendre (P1 (t)) which is defined as (t) · u (0). P1 (t) = u
(1)
Reorientation of end-to-end vectors could be one of the ways to characterize the PS whole chain movement and it is used here.
Fig. 2. Orientation correlation functions of PS end-to-end vectors from atomistic scale (left) and mesoscale (right) simulations in the 36PI and 24PS systems. At long times we can expect significant statistical error.
38
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
superatom chain into two parts, the middle and end part. The middle group includes the central three atoms or superatoms based on atomistic or mesoscale simulations, while the end group is composed by the three monomers closest to the ends on either side. The corresponding correlation functions are averaged to understand the heterogeneity along the chain in the dynamics of segments. In both atomistic and mesoscale models, the ends of the chains relax much faster than the middle ones as expected (Fig. 4). 3.3. Phase separation in PI–PS blends
Fig. 3. Rotation correlation function of PI (left) and PS (right) end-to-end vectors at various polymer systems.
Computer simulations can address not only the average properties of the system, but also the distribution over any observable of interest. The PS chains could, e.g., be classified as “fast”, “intermediate” and “slow” by splitting the total 24 chains into groups according to their rotation speed as seen from Fig. 2. We now perform a dynamic analysis for the atomistic simulations of PI–PS blends. We select systems of 58PI–2PS, 48PS, 36PI–24PS as examples. By comparing the average reorientations of PI and PS from 36PI–24PS, we confirmed the experimental suggestion [27] that the correlation times of PS are much slower than those of PI as shown in Fig. 3. The dynamics of polymer chains are not only determined by their own identity, but also strongly affected by their surrounding chains. By investigating the relaxation process of PI in 58PI–2PS and 36PI–24PS or PS in 48PS and 36PI–24PS, we validate that polymer chains movements do speed up or slow down depending on the average dynamics of neighboring chains and as shown above this leads to a wide variety of individual re-orientation functions. Fig. 3 clearly shows that the PI chains loose their initial orientation much faster with the increase of the PI concentration and PS chains become slower with the increase of PS concentration. The segmental motion of PS chains has also been characterized. We divide the atomistic 15 monomer or mesoscale 15
Using our mesoscale models we can now embark on a study of the phase behavior. We observe clear lamellae pattern for an equiconcentrated (by mass) mixture for 36PI–24PS30. The systems take on a cylindrical shape as the concentration of monomers changes to 60:40 and 66:34 in the cases of 36PI45–24PS30 and 36PI60–24PS30, respectively. The morphology changes into a spherical shape as the concentration goes to 72:28 as in the system of 36PI80–24PS30 [2,3]. To characterize the lamellar morphologies we calculate the density profiles of systems with a weight concentration of 50:50 in binary blends from both atomistic and mesoscale simulations. The chain lengths are 15, 30, and 60 monomers. We observe clearly phase separation behaviors in all cases, the only difference is the thickness of the interface region. As indicated in our previous work, a 15 monomer chain is too short to be classified as a polymer, it is only oligomers [1]. Polymers, at chain length of 30 or above, have to be used to investigate the characteristics of polymer blends. Fig. 5(b and c) are manifestations of a lamellar morphology that has formed during the phase separation as in two dimensions the densities are flat and the concentrations varies in one dimension. 3.4. Empirical mixing rules Let us now discuss the inapplicability of empirical mixing rules for this system. As mentioned above, we performed test simulations using algebraically and geometrically aver-
Fig. 4. Rotation correlation function of PS segmental vectors in 36PI15mer-24PS15mer from atomistic scale (left) and mesoscale (right).
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
39
Fig. 5. Density profile of 36PI–24PS at (a) chain length of 15 monomer from atomistic scale, (b) chain length of 30 monomer from mesoscale and (c) chain length of 60 monomer from mesoscale.
aged homo-interactions (PI–PI, PS–PS) as approximations for the hetero-interaction (PIPS), i.e., we approximated VPI–PS = 0.5(VPI–PI + VPSPS ) or VPI–PS = (VPI–PI + VPSPS )0.5 . At large distances where the VPS–PS is slightly negative we set VPI–PS = 0. Fig. 6 shows final results of systems which were set
up in morphologies equilibrated by the systematically mapped potentials. The empirical mixing rules clearly underestimate the phase separation tendency and lead to mixed systems in contrast to experimental results and the systematically mapped simulations. These comparison results clearly validate that an
Fig. 6. Snapshots of mixed systems using empirical mixing rules after being started in a morphology equilibrated by the systematic potential (a) 36PI6024PS60 using a geometrical averaging started from lamellae and (b) 36PI4524PS30 using an algebraic averaging started from a cylinder.
40
Q. Sun et al. / Fluid Phase Equilibria 261 (2007) 35–40
independent optimization of the hetero-interactions is inevitable in capturing the phase separation behaviors. 4. Conclusion Polymer melts and blends of chain lengths ranging from 7 up to 240 monomers with blends of various concentrations have been investigated from both atomistic and mesoscale models. We applied IBI in mapping the mesoscale models from atomistic simulations. With these mesoscale models, we reproduce the entanglement behaviors for longer PS chains and quantify the entanglement length of PS melt. Overall, the dynamics of PS chains is slower than that of PI and it speed up with the number of surrounding PS chains. There are clearly dynamic heterogeneity along the chain and among the chains. Phase separation of PI–PS blends is observed at chain length of 30 monomer and up. The morphology takes on a lamellar shape at concentration of 50:50, cylindrical at 60:40 and 66:34, or leads to a spherical shape at 72:28. The density profiles of blends are also characterized to provide quantitative results. It is noteworthy that in one of our earlier publications [2] we preliminarily suggested, based solely on RDFs, a cylindrical morphology for chains of length 80. After a more detailed analysis it turns out that that system is lamellar as well. It is clear that in principle the mesoscale potentials depend on concentration. Therefore, as we use the potential developed at a certain concentration we would expect that the boundaries between the different morphologies are only tentative. But it is clear that this approach is superior to the application of mixing rules as they lead to qualitatively wrong mixing behavior. Acknowledgments Financial support from the US Department of Energy, Office of Advanced Scientific Computing, under grant DE-FG0203ER25568 is gratefully acknowledged. References [1] Q. Sun, R. Faller, Crossover from unentangled to entangled dynamics in a systematically coarse-grained polystyrene melt, Macromolecules 39 (2) (2006) 812–820. [2] Q. Sun, R. Faller, Systematic coarse-graining of a polymer blend: polyisoprene and polystyrene, J. Chem. Theor. Comput. 2 (3) (2006) 607–615. [3] Q. Sun, R. Faller, Phase separation in polyisoprene/polystyrene blends by a systematically coarse-grained model, J. Chem. Phys. 126 (14) (2007) 144908. [4] W. Tsch¨op, K. Kremer, J. Batoulis, T. B¨urger, O. Hahn, Simulation of polymer melts. I. Coarse-graining for polycarbonates, Acta Polym. 49 (1998) 61–74. [5] J. Baschnagel, K. Binder, P. Doruker, A.A. Gusev, O. Hahn, K. Kremer, W.L. Mattice, F. M¨uller-Plathe, M. Murat, W. Paul, S. Santos, U.W. Suter, V. Tries, Bridging the gap between atomistic and coarse-grained models of polymers: Status and perspectives, Adv. Polym. Sci. 152 (2000) 141–156.
[6] H. Meyer, O. Biermann, R. Faller, D. Reith, F. M¨uller-Plathe, Coarse graining of nonbonded interparticle potentials using automatic simplex optimization to fit structural properties, J. Chem. Phys. 113 (15) (2000) 6264–6275. [7] R. Faller, F. M¨uller-Plathe, Multi-scale modeling of poly(isoprene) melts, Polymer 43 (2) (2002) 621–628. [8] R. Faller, Coarse-grain modeling of polymers, in: K. Lipkowitz, T. Cundari (Eds.), Reviews in Computational Chemistry, vol. 23, Wiley-VCH, 2007, pp. 233–262 (Chapter 4). [9] D. Reith, M. P¨utz, F. M¨uller-Plathe, Deriving effective mesoscale potentials from atomistic simulations, J. Comput. Chem. 24 (2003) 1624–1636. [10] G.S. Grest, M.-D. Lacasse, K. Kremer, A.M. Gupta, Efficient continuum model for simulating polymer blends and copolymers, J. Chem. Phys. 105 (23) (1996) 10583–10594. [11] K. Kremer, G.S. Grest, Dynamics of entangled linear polymer melts: a molecular-dynamics simulation, J. Chem. Phys. 92 (8) (1990) 5057–5086. [12] V.A. Harmandaris, V.G. Mavrantzas, D.N. Theodorou, M. Kr¨oger, J. ¨ Ramirez, H.C. Ottinger, D. Vlassopoulos, Crossover from the rouse to the entangled polymer melt regime: signals from long, detailed atomisitic molecular dynamics simulations, supported by rheological experiments, Macromolecules 36 (2003) 1376–1387. [13] S. Le´on, N. van der Vegt, L.D. Site, K. Kremer, Bisphenol a polycarbonate: entanglement analysis from coarse-grained md simulations, Macromolecules 38 (2005) 8078–8092. [14] Y. He, T.R. Lutz, M.D. Ediger, M. Pitsikalis, N. Hadjichristidis, E.D. von Meerwall, Miscible polyisoprene/polystyrene blends: distinct segmental dynamics but homogeneous terminal dynamics, Macromolecules 38 (2005) 6216–6226. [15] R. Koningsveld, W.J. MacKnight, Liquid–liquid phase separation in multicomponent polymer systems XXVII. Determination of the pair interaction function for polymer blends, Polym. Int. 44 (3) (1997) 356–364. [16] R. Faller, F. Muller-Plathe, M. Doxastakis, D. Theodorou, Local structure and dynamics of trans-polyisoprene oligomers, Macromolecules 34 (5) (2001) 1436–1448. [17] B. Hess, H. Bekker, H.J.C. Berendsen, J.G.E.M. Fraaije, LINCS: a linear constraint solver for molecular simulations, J. Comput. Chem. 18 (12) (1997) 1463–1472. [18] Q. Sun, R. Faller, Molecular dynamics of a polymer in mixed solvent: atactic polystyrene in a mixture of cyclohexane and n,n-dimethylformamide, J. Phys. Chem. B 109 (33) (2005) 15714–15723. [19] H.J.C. Berendsen, J.P.M. Postma, V. Gunsteren, D.W.F.A.J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (8) (1984) 3684–3690. [20] Q. Sun, R. Faller, Systematic coarse-graining of atomistic models for simulation of polymeric systems, Comput. Chem. Eng. 29 (11–12) (2005) 2380–2385. [21] R. Faller, D. Reith, Properties of polyisoprene model building in the melt and in solution, Macromolecules 36 (14) (2003) 5406–5414. [22] S. Hendrik Meyer, Private correspondence, Unpublished data, 2005. [23] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Oxford University Press Inc., 1987. [24] G.R. Strobl, The Physics of Polymers: Concepts for Understanding their Structures and Behavior, second ed., Springer, 1997. [25] L.J. Fetters, D.J. Lohse, D. Richter, T.A. Witten, A. Zirkel, Reviews connection between polymer molecular weight, density, chain dimensions, and melt viscoelastic properties, Macromolecules 27 (17) (1994) 4639–4647. [26] R. Faller, Correlation of static and dynamic inhomogeneities in polymer mixtures: a computer simulation of polyisoprene and polystyrene, Macromolecules 37 (3) (2004) 1095–1101. [27] Y. He, T.R. Lutz, M.D. Ediger, Segmental and terminal dynamics in miscible polymer blends, J. Chem. Phys. 119 (18) (2003) 9956–9965.