- Email: [email protected]
Primitive chain network model for block copolymers
Journal of Non-Crystalline Solids 352 (2006) 5001–5007 www.elsevier.com/locate/jnoncrysol
Primitive chain network model for block copolymers Yuichi Masubuchi a
a,b,*
, Giovanni Ianniruberto c, Francesco Greco d, Giuseppe Marrucci
c
Department of Organic and Polymer Materials Chemistry, Tokyo University of Agriculture and Technology, Tokyo 184-8588, Japan b Japan Science and Technology Agency, Tokyo 184-8588, Japan c Dipartimento di Ingegneria Chimica, Universita` degli studi di Napoli ‘‘Federico II’’, Piazzale Tecchio, 80-80125 Napoli, Italy d Instituto per i Materiali Compositi e Biomedici, CNR, Piazzale Tecchio, 80-80125 Napoli, Italy Available online 24 August 2006
Abstract A coarse-grained molecular simulation for block copolymers in the entangled state is proposed as an extension of the primitive chain network model. Polymers are represented as a sequence of segments between consecutive entanglements, the latter being modeled as sliplinks with other chains. Each sliplink connects two chains only, i.e., entanglements are taken as ‘binary’. The resulting 3D structure is a network of primitive chains, which makes our model different from other sliplink single-chain models, where the link to other chains is ‘virtual’. The 3D nature of our simulation makes it similar to, though considerably more coarse grained than, conventional molecular dynamics simulations. Because of the 3D space assignment of the polymers, monomeric density fields can be defined, and interactions due to different chemistry of the monomers can be accounted for, similarly to density field calculations. Polymer motion is described both by the 3D motion of sliplinks, and by the 1D transport of monomers along the primitive chain, while network topological rearrangement occurs due to chain-end hooking and unhooking processes. Each kinetic equation accounts for elastic forces along the chains, field forces arising from density gradients, and thermal random forces. In this paper, the primitive chain network model was modified (i) in the procedure of network rearrangement to account for the different chemistries in the copolymer, and (ii) in some details of the kinetic equations whenever the boundary between blocks is involved. We report results for diblock copolymers where for simplicity all relevant properties of the two monomers are the same, except for the interactions. Simulations reasonably reproduce the micro-phase formation process and the phase diagram for well entangled copolymers with a low calculation cost. 2006 Elsevier B.V. All rights reserved. PACS: 83.10.Kn; 83.10.Mj; 82.35.Jk Keywords: Modeling and simulation; Polymers and organics; Rheology
1. Introduction Because of the increasing importance of block copolymers [1], various simulation methods have been proposed to predict microphase structures and polymer dynamics. For the equilibrium microphase structure, density field based methods [2] and lattice Monte Carlo simulations [3,4] have been reported. For polymer dynamics and rheo* Corresponding author. Address: Department of Organic and Polymer Materials Chemistry, Tokyo University of Agriculture and Technology, Tokyo 184-8588, Japan. Tel./fax: +81 42 388 7058. E-mail address: [email protected] (Y. Masubuchi).
0022-3093/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2006.01.144
logical properties, microscopic [5] and coarse-grained [6] molecular dynamics and dissipative particle dynamics [7] have been used. In some cases, a combination of different methods is preferred. For example, Aoyagi et al. [8] used a self-consistent field calculation to obtain the equilibrium structures, then performing coarse-grained molecular dynamics on those structures to calculate dynamics and rheology. The above mentioned methods, certainly useful for prediction of equilibrium structures, can in fact afford the dynamics only for the case of short chains or for the short time response. Block copolymers in the entangled state remain a challenging problem both because of the specific
5002
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
features of the entangled polymer dynamics and for the extremely long relaxation time. For example, the coarsegrained molecular dynamics is not practically used for entangled systems, since few entanglements per chain is an effective upper limit due to computational time limitations. Although dissipative particle methods are more highly coarse-grained, they remain restricted to unentangled situations due to the truncated interaction potentials. To treat the entangled polymer dynamics, sliplink models [9–13] have been developed based on the established reptation concept [14]. Sliplink models could predict rheology and polymer dynamics up to the longest relaxation time for entangled homopolymer systems, and also an extended model for polymer blends has been proposed [15]. However, no sliplink model has so far been used for block copolymers. A promising option is a combination of reptation and density field ideas. Recently, Shima et al. [16] proposed a model incorporating reptation dynamics into the monomer flux of density field based methods. Taking into account biased reptation, chain retraction, and constraint release, they showed in particular that the specific features of entangled brushes were qualitatively reproduced. Although probably obtainable in some similar way, results for entangled block copolymers from the method of Shima et al. [16] have not so far been reported. In this work, we extend the primitive chain network model [10,15] to the case of block copolymers. In the primitive chain network model, differently from other sliplink models, the chains are assigned in real space and distributed in a simulation box, similarly to conventional molecular dynamics simulations. Due to the real space assignment, force balance on entanglements and (more importantly for the case considered here) field forces due to density gradients can be considered, in addition to the known mechanisms of entangled polymer dynamics, i.e., reptation and contour length fluctuation [14], thermal [17] and convective [18] constraint release. In our previous paper [15], we extended the model proposed in [10] to polymer blends by introducing interactions between different chemicals through the field force. In the present work on block copolymers, we further modify the model by accounting for the fact that the chain-end hooking process that creates new entanglements should be biased by the interaction. We also modify the expression for the tension in the mixed-chemistry chain segment that includes the boundary between the blocks.
links connecting chains in pair, so that the chains (though keeping their internal connectivity) effectively form an infinite network of primitive paths, with a mesh size equal to the mean segment length. Fig. 1 schematically represents the variables of state describing position R of sliplinks and number of monomers in chain segments between sliplinks. Kinetic equations for these variables will be presented below. Although chains slide through them, sliplinks are permanent, except at chain ends, where they can be abandoned and therefore vanish. The model specifies how existing sliplinks disappear and new ones are created at chain ends. All together, the dynamics of the chains is described by (i) motion of the sliplinks in 3D space (similarly to the crosslinks of a rubber network), (ii) monomer transport along the primitive path, hence through the sliplinks (reptation plus fluctuations plus deformation-induced retraction), and (iii) network rearrangement due to sliplink destruction and reformation at chain ends (constraint renewal). It is worth emphasizing that, consistently with the level of coarse-graining adopted in our model, chain segments between sliplinks are taken to be phantom and hence, differently from more detailed entanglement models [5,19], two-body interactions are not explicitly accounted for. Implicitly, however, interactions are somehow included in our model through a mean-field density-preserving potential. 2.1. Sliplink motion The equation of motion of the sliplinks is written as: ðfa þ fb ÞðR_ k RÞ " # 1 ri ri1 1 rj rj1 ¼ 3kT 2 þ 2 ba ni ni1 bb nj nj1 rðn0a la þ n0b lb Þ þ F:
ð1Þ
2. Model The model is essentially the same of the primitive chain network model for blends [15], with the improvements anticipated in the Introduction and detailed later. In the primitive chain network model, polymer chains are represented by connected segments, each between consecutive entanglements, and are dispersed in a simulation box. Entanglements among polymers are represented as slip-
Fig. 1. Schematic representation of the sliplink model with relevant state variables. Detailed description is given in the text.
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
The left hand side of this equation is the friction force on the sliplink, where the dot on sliplink position R indicates time derivative (i.e., sliplink velocity), k is the velocity gradient tensor of the continuum, and f is the friction coefficient of the pair of segments linked by the sliplink. Suffices a and b indicate different types of monomers, but of course we can have a = b if the sliplink connects two segments with the same chemistry (compare Fig. 1). The first term on the right hand side of Eq. (1) is the elastic contribution of the segments pulling on the sliplink. Here, kT is thermal energy, b is monomer length, and r and n are chain-segment end-to-end vector and monomer number, respectively. Indices on r and n specify the position of a segment along its own chain, and vector r is oriented along increasing values of the index. If segment i happens to be of mixed chemistry (i.e., if it includes the block boundary), Gaussian statistics gives the elastic force Fel i in such a segment as: Fel i ¼
3kT ri : þ b2b nib
ð2Þ
b2a nia
Note that, although the simulation does not monitor the position in space of the block boundary, the number of monomers for each chemical species present in the segment is monitored, the block boundary sliding through the sliplinks together with the chain as detailed below. Note further that also the friction coefficient of the mixed segment can account for the different chemistry, linearly with the monomer number in the segment. The second term on the right hand side of Eq. (1) is the field force due to the gradient of the monomer chemical potential l, itself derived from the free energy: A ¼ Amix þ Avol ;
ð3Þ
where Amix is the enthalpic contribution of the Flory–Huggins expression for binary polymer blends [20], and Avol is a phenomenological free energy that accounts, as previously mentioned, for the constant density constraint: Amix =kT ¼ v/a /b ; 8 2 < / e h/i 1 Avol =kT ¼ : 0
ð4Þ for / > h/i;
ð5Þ
for / 6 h/i:
Here / is the monomer volume fraction, v is the Flory interaction parameter, and e is the phenomenological modulus for incompressibility [15]. We recall that, since our chain segments are phantom, without Avol chains and sliplinks tend to cluster during the simulation run. In the gradient term of Eq. (1) the monomer chemical potential is multiplied to the monomer number in the chain segment, for simplicity taken as the equilibrium mean value n0. Finally, the third term of Eq. (1) is a Gaussian random force obeying fluctuation–dissipation theorem [14]: hFi ¼ 0; Fi ðtÞ Fj ðt0 Þ ¼ 6kT ðfa þ fb Þdij dðt t0 Þ: ð6Þ
5003
2.2. Monomer transport along the primitive path Monomer sliding through the sliplinks is described by the following kinetic equation for the rate of change of monomer number n_ in the two chain segments across a sliplink (n_ only differing in sign in the two of them): n_ 3kT ri ri1 ð7Þ n0a rla þ f : fa ¼ 2 q ba ni ni1 The left hand side of this equation is the scalar friction force associated to sliding, q being the average linear monomer density in the two segments involved, i.e., the mean local density along the contour of the primitive path: 1 n n q¼ þ : ð8Þ 2 r i r i1 The first term on the right hand side of Eq. (7) is the elastic tension difference (here r is the magnitude of r). Consistently with the above, if one of the two segments is of mixed chemistry, the elastic tension is calculated through Eq. (2). The second term on the right hand side is the field term derived from the free energy given by Eq. (3), where now the chemical potential gradient is calculated along the primitive path rather than in 3D space. Finally the third term is the 1D Gaussian force that obeys hf i ¼ 0;
hfi ðtÞfj ðt0 Þi ¼ 2kT fa dij dðt t0 Þ:
ð9Þ
2.3. Network topology rearrangement As mentioned above, chain ends govern creation and destruction of sliplinks. The rule for topological changes is based upon monomer number in the end segments of the primitive chains. Specifically, if the monomer number in a chain end falls below a certain minimum value, the end sliplink is destroyed (since we assume that the chain end instantly slips through it). As a consequence, the partner chain in the 3D network also looses that constraint (constraint release). Conversely, if the monomer number in a chain end exceeds a maximum value, a new sliplink is created on the end segment by hooking one of surrounding segments. The number window employed is 0:5n0 < n < 1:5n0 :
ð10Þ
In our previous paper on blends [15], sliplink creation was performed by randomly hooking surrounding segments independently of their chemistry. However, it appears reasonable to assume that interactions also come into play in the selection of the partner in the hooking process, since contacts with like segments are systematically preferred over those with different ones. To introduce the effect of the interaction between different chemicals in the hooking process, a sort of Monte Carlo procedure is employed for the selection of the partner [21]. When the randomly selected segment is of the same chemistry of the hooking end, the new entanglement is accepted as in the conventional procedure. Conversely, if
5004
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
the candidate partner has a different chemistry, acceptance of the hooking is subjected to the Boltzmann factor: p ¼ expðvn0 Þ;
ð11Þ
where n0 is some average value between n0a and n0b . If the hooking is rejected, the selection procedure of the partner segment is repeated until, one way or the other, hooking is accepted. 3. Simulation details In this work we have examined entangled linear diblock copolymers. Since we are here only interested in qualitative aspects, we have assumed the simplest possible situation, i.e., that the two chemicals have the same monomer length (actually Kuhn segment length) b, the same average monomer number n0 per primitive chain p segment, and hence the same entanglement spacing a = b n0. We also assume equal friction coefficient f, and equal monomeric volume. Calculations were performed by choosing a as the unit length, kT as the unit energy, and a2f/6kT as the unit time. Also, monomer numbers are normalized to n0, i.e., in the calculations n0 is taken as unity. Real units will obviously depend on the specific polymer. By way of example, for a polystyrene melt at 160 C, the number of monomers n0 between consecutive entanglements in our code corresponds to a molecular mass of ca. 8300 g/mol, from which ˚ is obtained, while the unit time a value for a of ca. 60 A comes out equal to ca. 0.01 s [22]. Euler numerical integration with a (non-dimensional) time step of Dt = 0.01 was adopted. A periodic box with a size of either 123 or 163 was used, depending on the polymer length. The number density of primitive chain segments was fixed at 10, a typical value for melts. To calculate the field term, the local chemical potential was obtained by counting monomers within subdivisions of the simulation box. To this purpose we used cubes with size of 13. Finally, the phenomenological incompressibility parameter e in Fvol was set to 0.5. 4. Results Fig. 2 shows typical snapshots of the simulation box for diblock copolymers with a total monomer number N = 40n0, i.e., with ca. 40 entanglements per chain. The block ratio, coinciding in our case with the volume fraction, is /a = /b = 0.5, and the interaction parameter is vn0 = 0.5. The figure shows the situation reached in 100 000 time units after quenching from the homogeneous state. Fig. 2(a) shows one of the molecules in the system to compare chain size to that of the 163 box. Blue and yellow colors of the segments indicate the two different chemistries in the test chain.1 Red thin lines show segments of 1 For interpretation of color in Fig. 2, the reader is referred to the web version of this article.
Fig. 2. Snapshots of an entangled (N = 40n0) diblock copolymer simulation with a block ratio of 0.5 and vn0 = 0.5, taken at 100 000 time units after quenching from the homogeneous state. (a) A typical molecule of the system, (b) sliplink distribution in the simulation box, and (c) density field.
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
surrounding chains entangled with the test chain. Fig. 2(b) shows the distribution in space of the sliplinks in the whole box. Yellow and blue cubes correspond to sliplinks linking segments with the same chemistry. Conversely, green cubes indicate sliplinks joining chemically different segments. A lamellar-like structure can be recognized, as one would
5005
expect for the equi-diblock copolymer. Green cubes are seen to concentrate in the phase boundary region. Finally Fig. 2(c) shows with a color scale the monomer density in the mesh used in the chemical potential gradient calculations, so as to appreciate the spatial resolution of the field force.
Fig. 3. Time evolution of the system shown in Fig. 2. Time units after quenching are: (a) 0; (b) 2000; (c) 4000; (d) 10 000; (e) 20 000; and (f) 40 000. The 100 000 snapshot is in Fig. 2(b).
5006
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
Fig. 3 shows the evolution in time of the system presented in Fig. 2. In the early stages of the phase separation (see Fig. 3(a)–(c)), polymer reptation does not occur over significant distances, and hence formation of the microphase structure occurs through local motions of both chains and sliplinks. In fact, the longest relaxation time of a polymer having N = 40n0 is around 5000. On the other hand, in later stages (see Fig. 3(d)–(f)) the liquid nature of the system emerges, allowing for diffusion over long distances. Motion and deformation of the whole microphase structure can then be observed. Fig. 4 shows the phase diagram for both N = 40n0 and N = 20n0 diblock copolymers. Calculations were performed with a 163 simulation box for N = 40n0, and with a 123 box for N = 20n0. Predictions of a self-consistent field calculation (SCF) [23] are also presented for comparison. The phase boundary between the disordered phase and the ordered ones is reasonably reproduced by our simulations. It must be said, however, that for the case N = 10n0 no ordered phase could be observed. This is most probably due to the poor resolution of the mesh used for chemical potential gradient calculations, which for a 10n0 chain is indeed comparable with the block size. On the other hand, finer meshes are not practicable because the statistics (with the standard segment density of 10) are not sufficient. It so appears that the results obtained so far are valid for well entangled polymers only, which is indeed the main aim of the present contribution. Finally, Fig. 5 shows the kind of microphase structures obtained, depending on block ratio. Reasonable reproductions of the sphere and cylinder phases can be recognized. However, in the present work the size of the simulation box was fixed independently of the equilibrium microphase structure to be eventually reached, and was not optimized.
80 SCF N = 40n0 N = 20n0
60
χN
Hence the ordered phases, as obtained in classical density field calculations, were not exactly reproduced. 5. Discussion
S
C
L
40 G 20
0
Fig. 5. Snapshots after 100 000 time units for diblock copolymers with the same molecular mass as that in Fig. 1 (N = 40n0), but with different block ratios and interaction parameters: (a) / = 0.1 and vn0 = 2.0; (b) / = 0.3 and vn0 = 2.0.
DISORDERED 0.0
0.1
0.2 0.3 block ratio
0.4
0.5
Fig. 4. The phase diagram for two diblock copolymers with N = 40n0 and N = 20n0. Filled symbols indicate microphase separation and open ones correspond to the disordered state. SCF results [23] are also presented for comparison (lines) with the indication of the ordered phase structure: (S) sphere, (C) cylinder, (L) lamella, and (G) gyroid.
In this section, we discuss how the proposed model compares with, and possibly complement, the conventional methods employed for block copolymers. It is emphasized that our method is the single molecularly-based model able to describe with acceptable computational cost the global motion of well-entangled block copolymers. Indeed computations were performed with a Linux workstation equipped with a 3.0 GHz Pentium 4 processor with 1 GB main memory. With that computer, and for the 163 simulation box filled with the typical segment density of 10 (i.e., for a total number of segments of 10 · 163 = 4.1 · 104), ca. 16 h are required for a run of 5000 time units, roughly corresponding to the longest relaxation time of chains with 40 entanglements. It is also noteworthy that the calculation time is proportional to
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
the system size since the model does not account for twobody non-bond interactions, which are typically the most expensive in conventional molecular dynamics as their number increases with the square of the system size. However, as mentioned before, our computations have shown some difficulties with marginally entangled chains because of the resolution of the chemical potential field. Also in view of the fact that our model certainly does not apply to unentangled chains, multi-scale simulations are required if one needs to explore the whole range of molecular weights. Of course, our computations of the dynamical evolution are not fast if compared with the density field calculations used to obtain equilibrium structures. For instance, a recently reported simulation [24] takes 3 h (with a similar computer to that used in the present work) to obtain the equilibrium structure of a mixture of block-copolymers and linear polymers in a box of 643 lattice points. On the other hand, since most density-field-based dynamic simulations do not account for entanglements, a suitable combination of different methods complementing each other seems promising. Indeed, the equilibrium structure quickly obtained with the density field theories certainly remains useful because the present method cannot easily reproduce the exact ordered structures of the microphases. 6. Conclusions A simulation method for entangled block copolymers was here proposed as an extension of the primitive chain network model [10,15]. Predictions for the phase diagram are in qualitative agreement with the conventional methods, but phase separation kinetics is also described. Hence, the present method may occupy an important niche among the existing ones, falling somewhere in between molecular dynamics [6], density field calculations [2], and singlechain sliplink simulations [9,11]. Detailed assessment of the method by quantitative comparison with experiments on rheology, phase separation kinetics, molecular confor-
5007
mation and molecular motion, etc. is in progress, and will be presented elsewhere. Acknowledgement This work was supported in part by a Grant in Aid for Scientific Research [KAKENHI] 2005 No. 17750200 from the Ministry of Education, Culture, Sports, Science and Technology of Japan. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
G. Srinivas, D.E. Discher, M.L. Klein, Nature Mater. 3 (2005) 638. G.H. Fredrickson, V. Ganesan, F. Drolet, Macromolecules 35 (2002) 16. K. Binder, Comp. Phys. Commun. 147 (2002) 22. T. Dotera, A. Hatano, J. Chem. Phys. 105 (1996) 8413. A. Terzis, D.N. Theodorou, A. Stroeks, Macromolecules 35 (2002) 508. M. Murat, G.S. Grest, K. Kremer, Macromolecules 32 (1999) 595. P. Ahlrichs, B. Dunweg, J. Chem. Phys. 111 (1999) 8225. T. Aoyagi, T. Honda, M. Doi, J. Chem. Phys. 117 (2002) 8153. C.C. Hua, J. Schieber, J. Chem. Phys. 109 (2001) 10018. Y. Masubuchi, J.-I. Takimoto, K. Koyama, G. Ianniruberto, G. Marrucci, F. Greco, J. Chem. Phys. 115 (2001) 4387. J.-I. Takimoto, M. Doi, Philos. Trans. Royal Soc. London A 361 (2003) 641. J.D. Schieber, J. Neergaard, S. Gupta, J. Rheol. 47 (2003) 213. A.E. Likhtman, Macromolecules 38 (2005) 6128. M. Doi, S.F. Edwards, The Theory of Polymer Dynamics, Clarendon, Oxford, 1986. Y. Masubuchi, G. Ianniruberto, F. Greco, G. Marrucci, Model. Simul. Mater. Sci. Eng. 12 (2004) S91. T. Shima, H. Kuni, Y. Okabe, M. Doi, X.F. Yuan, T. Kawakatsu, Macromolecules 36 (2003) 9199. J. des Cloizeaux, Europhys. Lett. 5 (1988) 437. G. Marrucci, J. Non-Newtonian Fluid Mech. 62 (1996) 279. J.T. Padding, W.J. Briels, J. Chem. Phys. 117 (2002) 925. P.J. Flory, Principles of Polymer Chemistry, Cornell University Press, Ithaca, 1953. Y. Masubuchi, G. Ianniruberto, F. Greco, G. Marrucci, KoreaAustralia Rheology J. 18 (2006) 99. T. Yaoita, T. Isaki, Y. Masubuchi, H. Watanabe, G. Ianniruberto, F. Greco, G. Marrucci, J. Chem. Phys. 121 (2004) 12650. M.W. Matsen, F.S. Bates, Macromolecules 29 (1996) 1091. T. Uneyama, M. Doi, Macromolecules 38 (2005) 196.
Primitive chain network model for block copolymers Yuichi Masubuchi a
a,b,*
, Giovanni Ianniruberto c, Francesco Greco d, Giuseppe Marrucci
c
Department of Organic and Polymer Materials Chemistry, Tokyo University of Agriculture and Technology, Tokyo 184-8588, Japan b Japan Science and Technology Agency, Tokyo 184-8588, Japan c Dipartimento di Ingegneria Chimica, Universita` degli studi di Napoli ‘‘Federico II’’, Piazzale Tecchio, 80-80125 Napoli, Italy d Instituto per i Materiali Compositi e Biomedici, CNR, Piazzale Tecchio, 80-80125 Napoli, Italy Available online 24 August 2006
Abstract A coarse-grained molecular simulation for block copolymers in the entangled state is proposed as an extension of the primitive chain network model. Polymers are represented as a sequence of segments between consecutive entanglements, the latter being modeled as sliplinks with other chains. Each sliplink connects two chains only, i.e., entanglements are taken as ‘binary’. The resulting 3D structure is a network of primitive chains, which makes our model different from other sliplink single-chain models, where the link to other chains is ‘virtual’. The 3D nature of our simulation makes it similar to, though considerably more coarse grained than, conventional molecular dynamics simulations. Because of the 3D space assignment of the polymers, monomeric density fields can be defined, and interactions due to different chemistry of the monomers can be accounted for, similarly to density field calculations. Polymer motion is described both by the 3D motion of sliplinks, and by the 1D transport of monomers along the primitive chain, while network topological rearrangement occurs due to chain-end hooking and unhooking processes. Each kinetic equation accounts for elastic forces along the chains, field forces arising from density gradients, and thermal random forces. In this paper, the primitive chain network model was modified (i) in the procedure of network rearrangement to account for the different chemistries in the copolymer, and (ii) in some details of the kinetic equations whenever the boundary between blocks is involved. We report results for diblock copolymers where for simplicity all relevant properties of the two monomers are the same, except for the interactions. Simulations reasonably reproduce the micro-phase formation process and the phase diagram for well entangled copolymers with a low calculation cost. 2006 Elsevier B.V. All rights reserved. PACS: 83.10.Kn; 83.10.Mj; 82.35.Jk Keywords: Modeling and simulation; Polymers and organics; Rheology
1. Introduction Because of the increasing importance of block copolymers [1], various simulation methods have been proposed to predict microphase structures and polymer dynamics. For the equilibrium microphase structure, density field based methods [2] and lattice Monte Carlo simulations [3,4] have been reported. For polymer dynamics and rheo* Corresponding author. Address: Department of Organic and Polymer Materials Chemistry, Tokyo University of Agriculture and Technology, Tokyo 184-8588, Japan. Tel./fax: +81 42 388 7058. E-mail address: [email protected] (Y. Masubuchi).
0022-3093/$ - see front matter 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.jnoncrysol.2006.01.144
logical properties, microscopic [5] and coarse-grained [6] molecular dynamics and dissipative particle dynamics [7] have been used. In some cases, a combination of different methods is preferred. For example, Aoyagi et al. [8] used a self-consistent field calculation to obtain the equilibrium structures, then performing coarse-grained molecular dynamics on those structures to calculate dynamics and rheology. The above mentioned methods, certainly useful for prediction of equilibrium structures, can in fact afford the dynamics only for the case of short chains or for the short time response. Block copolymers in the entangled state remain a challenging problem both because of the specific
5002
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
features of the entangled polymer dynamics and for the extremely long relaxation time. For example, the coarsegrained molecular dynamics is not practically used for entangled systems, since few entanglements per chain is an effective upper limit due to computational time limitations. Although dissipative particle methods are more highly coarse-grained, they remain restricted to unentangled situations due to the truncated interaction potentials. To treat the entangled polymer dynamics, sliplink models [9–13] have been developed based on the established reptation concept [14]. Sliplink models could predict rheology and polymer dynamics up to the longest relaxation time for entangled homopolymer systems, and also an extended model for polymer blends has been proposed [15]. However, no sliplink model has so far been used for block copolymers. A promising option is a combination of reptation and density field ideas. Recently, Shima et al. [16] proposed a model incorporating reptation dynamics into the monomer flux of density field based methods. Taking into account biased reptation, chain retraction, and constraint release, they showed in particular that the specific features of entangled brushes were qualitatively reproduced. Although probably obtainable in some similar way, results for entangled block copolymers from the method of Shima et al. [16] have not so far been reported. In this work, we extend the primitive chain network model [10,15] to the case of block copolymers. In the primitive chain network model, differently from other sliplink models, the chains are assigned in real space and distributed in a simulation box, similarly to conventional molecular dynamics simulations. Due to the real space assignment, force balance on entanglements and (more importantly for the case considered here) field forces due to density gradients can be considered, in addition to the known mechanisms of entangled polymer dynamics, i.e., reptation and contour length fluctuation [14], thermal [17] and convective [18] constraint release. In our previous paper [15], we extended the model proposed in [10] to polymer blends by introducing interactions between different chemicals through the field force. In the present work on block copolymers, we further modify the model by accounting for the fact that the chain-end hooking process that creates new entanglements should be biased by the interaction. We also modify the expression for the tension in the mixed-chemistry chain segment that includes the boundary between the blocks.
links connecting chains in pair, so that the chains (though keeping their internal connectivity) effectively form an infinite network of primitive paths, with a mesh size equal to the mean segment length. Fig. 1 schematically represents the variables of state describing position R of sliplinks and number of monomers in chain segments between sliplinks. Kinetic equations for these variables will be presented below. Although chains slide through them, sliplinks are permanent, except at chain ends, where they can be abandoned and therefore vanish. The model specifies how existing sliplinks disappear and new ones are created at chain ends. All together, the dynamics of the chains is described by (i) motion of the sliplinks in 3D space (similarly to the crosslinks of a rubber network), (ii) monomer transport along the primitive path, hence through the sliplinks (reptation plus fluctuations plus deformation-induced retraction), and (iii) network rearrangement due to sliplink destruction and reformation at chain ends (constraint renewal). It is worth emphasizing that, consistently with the level of coarse-graining adopted in our model, chain segments between sliplinks are taken to be phantom and hence, differently from more detailed entanglement models [5,19], two-body interactions are not explicitly accounted for. Implicitly, however, interactions are somehow included in our model through a mean-field density-preserving potential. 2.1. Sliplink motion The equation of motion of the sliplinks is written as: ðfa þ fb ÞðR_ k RÞ " # 1 ri ri1 1 rj rj1 ¼ 3kT 2 þ 2 ba ni ni1 bb nj nj1 rðn0a la þ n0b lb Þ þ F:
ð1Þ
2. Model The model is essentially the same of the primitive chain network model for blends [15], with the improvements anticipated in the Introduction and detailed later. In the primitive chain network model, polymer chains are represented by connected segments, each between consecutive entanglements, and are dispersed in a simulation box. Entanglements among polymers are represented as slip-
Fig. 1. Schematic representation of the sliplink model with relevant state variables. Detailed description is given in the text.
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
The left hand side of this equation is the friction force on the sliplink, where the dot on sliplink position R indicates time derivative (i.e., sliplink velocity), k is the velocity gradient tensor of the continuum, and f is the friction coefficient of the pair of segments linked by the sliplink. Suffices a and b indicate different types of monomers, but of course we can have a = b if the sliplink connects two segments with the same chemistry (compare Fig. 1). The first term on the right hand side of Eq. (1) is the elastic contribution of the segments pulling on the sliplink. Here, kT is thermal energy, b is monomer length, and r and n are chain-segment end-to-end vector and monomer number, respectively. Indices on r and n specify the position of a segment along its own chain, and vector r is oriented along increasing values of the index. If segment i happens to be of mixed chemistry (i.e., if it includes the block boundary), Gaussian statistics gives the elastic force Fel i in such a segment as: Fel i ¼
3kT ri : þ b2b nib
ð2Þ
b2a nia
Note that, although the simulation does not monitor the position in space of the block boundary, the number of monomers for each chemical species present in the segment is monitored, the block boundary sliding through the sliplinks together with the chain as detailed below. Note further that also the friction coefficient of the mixed segment can account for the different chemistry, linearly with the monomer number in the segment. The second term on the right hand side of Eq. (1) is the field force due to the gradient of the monomer chemical potential l, itself derived from the free energy: A ¼ Amix þ Avol ;
ð3Þ
where Amix is the enthalpic contribution of the Flory–Huggins expression for binary polymer blends [20], and Avol is a phenomenological free energy that accounts, as previously mentioned, for the constant density constraint: Amix =kT ¼ v/a /b ; 8 2 < / e h/i 1 Avol =kT ¼ : 0
ð4Þ for / > h/i;
ð5Þ
for / 6 h/i:
Here / is the monomer volume fraction, v is the Flory interaction parameter, and e is the phenomenological modulus for incompressibility [15]. We recall that, since our chain segments are phantom, without Avol chains and sliplinks tend to cluster during the simulation run. In the gradient term of Eq. (1) the monomer chemical potential is multiplied to the monomer number in the chain segment, for simplicity taken as the equilibrium mean value n0. Finally, the third term of Eq. (1) is a Gaussian random force obeying fluctuation–dissipation theorem [14]: hFi ¼ 0; Fi ðtÞ Fj ðt0 Þ ¼ 6kT ðfa þ fb Þdij dðt t0 Þ: ð6Þ
5003
2.2. Monomer transport along the primitive path Monomer sliding through the sliplinks is described by the following kinetic equation for the rate of change of monomer number n_ in the two chain segments across a sliplink (n_ only differing in sign in the two of them): n_ 3kT ri ri1 ð7Þ n0a rla þ f : fa ¼ 2 q ba ni ni1 The left hand side of this equation is the scalar friction force associated to sliding, q being the average linear monomer density in the two segments involved, i.e., the mean local density along the contour of the primitive path: 1 n n q¼ þ : ð8Þ 2 r i r i1 The first term on the right hand side of Eq. (7) is the elastic tension difference (here r is the magnitude of r). Consistently with the above, if one of the two segments is of mixed chemistry, the elastic tension is calculated through Eq. (2). The second term on the right hand side is the field term derived from the free energy given by Eq. (3), where now the chemical potential gradient is calculated along the primitive path rather than in 3D space. Finally the third term is the 1D Gaussian force that obeys hf i ¼ 0;
hfi ðtÞfj ðt0 Þi ¼ 2kT fa dij dðt t0 Þ:
ð9Þ
2.3. Network topology rearrangement As mentioned above, chain ends govern creation and destruction of sliplinks. The rule for topological changes is based upon monomer number in the end segments of the primitive chains. Specifically, if the monomer number in a chain end falls below a certain minimum value, the end sliplink is destroyed (since we assume that the chain end instantly slips through it). As a consequence, the partner chain in the 3D network also looses that constraint (constraint release). Conversely, if the monomer number in a chain end exceeds a maximum value, a new sliplink is created on the end segment by hooking one of surrounding segments. The number window employed is 0:5n0 < n < 1:5n0 :
ð10Þ
In our previous paper on blends [15], sliplink creation was performed by randomly hooking surrounding segments independently of their chemistry. However, it appears reasonable to assume that interactions also come into play in the selection of the partner in the hooking process, since contacts with like segments are systematically preferred over those with different ones. To introduce the effect of the interaction between different chemicals in the hooking process, a sort of Monte Carlo procedure is employed for the selection of the partner [21]. When the randomly selected segment is of the same chemistry of the hooking end, the new entanglement is accepted as in the conventional procedure. Conversely, if
5004
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
the candidate partner has a different chemistry, acceptance of the hooking is subjected to the Boltzmann factor: p ¼ expðvn0 Þ;
ð11Þ
where n0 is some average value between n0a and n0b . If the hooking is rejected, the selection procedure of the partner segment is repeated until, one way or the other, hooking is accepted. 3. Simulation details In this work we have examined entangled linear diblock copolymers. Since we are here only interested in qualitative aspects, we have assumed the simplest possible situation, i.e., that the two chemicals have the same monomer length (actually Kuhn segment length) b, the same average monomer number n0 per primitive chain p segment, and hence the same entanglement spacing a = b n0. We also assume equal friction coefficient f, and equal monomeric volume. Calculations were performed by choosing a as the unit length, kT as the unit energy, and a2f/6kT as the unit time. Also, monomer numbers are normalized to n0, i.e., in the calculations n0 is taken as unity. Real units will obviously depend on the specific polymer. By way of example, for a polystyrene melt at 160 C, the number of monomers n0 between consecutive entanglements in our code corresponds to a molecular mass of ca. 8300 g/mol, from which ˚ is obtained, while the unit time a value for a of ca. 60 A comes out equal to ca. 0.01 s [22]. Euler numerical integration with a (non-dimensional) time step of Dt = 0.01 was adopted. A periodic box with a size of either 123 or 163 was used, depending on the polymer length. The number density of primitive chain segments was fixed at 10, a typical value for melts. To calculate the field term, the local chemical potential was obtained by counting monomers within subdivisions of the simulation box. To this purpose we used cubes with size of 13. Finally, the phenomenological incompressibility parameter e in Fvol was set to 0.5. 4. Results Fig. 2 shows typical snapshots of the simulation box for diblock copolymers with a total monomer number N = 40n0, i.e., with ca. 40 entanglements per chain. The block ratio, coinciding in our case with the volume fraction, is /a = /b = 0.5, and the interaction parameter is vn0 = 0.5. The figure shows the situation reached in 100 000 time units after quenching from the homogeneous state. Fig. 2(a) shows one of the molecules in the system to compare chain size to that of the 163 box. Blue and yellow colors of the segments indicate the two different chemistries in the test chain.1 Red thin lines show segments of 1 For interpretation of color in Fig. 2, the reader is referred to the web version of this article.
Fig. 2. Snapshots of an entangled (N = 40n0) diblock copolymer simulation with a block ratio of 0.5 and vn0 = 0.5, taken at 100 000 time units after quenching from the homogeneous state. (a) A typical molecule of the system, (b) sliplink distribution in the simulation box, and (c) density field.
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
surrounding chains entangled with the test chain. Fig. 2(b) shows the distribution in space of the sliplinks in the whole box. Yellow and blue cubes correspond to sliplinks linking segments with the same chemistry. Conversely, green cubes indicate sliplinks joining chemically different segments. A lamellar-like structure can be recognized, as one would
5005
expect for the equi-diblock copolymer. Green cubes are seen to concentrate in the phase boundary region. Finally Fig. 2(c) shows with a color scale the monomer density in the mesh used in the chemical potential gradient calculations, so as to appreciate the spatial resolution of the field force.
Fig. 3. Time evolution of the system shown in Fig. 2. Time units after quenching are: (a) 0; (b) 2000; (c) 4000; (d) 10 000; (e) 20 000; and (f) 40 000. The 100 000 snapshot is in Fig. 2(b).
5006
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
Fig. 3 shows the evolution in time of the system presented in Fig. 2. In the early stages of the phase separation (see Fig. 3(a)–(c)), polymer reptation does not occur over significant distances, and hence formation of the microphase structure occurs through local motions of both chains and sliplinks. In fact, the longest relaxation time of a polymer having N = 40n0 is around 5000. On the other hand, in later stages (see Fig. 3(d)–(f)) the liquid nature of the system emerges, allowing for diffusion over long distances. Motion and deformation of the whole microphase structure can then be observed. Fig. 4 shows the phase diagram for both N = 40n0 and N = 20n0 diblock copolymers. Calculations were performed with a 163 simulation box for N = 40n0, and with a 123 box for N = 20n0. Predictions of a self-consistent field calculation (SCF) [23] are also presented for comparison. The phase boundary between the disordered phase and the ordered ones is reasonably reproduced by our simulations. It must be said, however, that for the case N = 10n0 no ordered phase could be observed. This is most probably due to the poor resolution of the mesh used for chemical potential gradient calculations, which for a 10n0 chain is indeed comparable with the block size. On the other hand, finer meshes are not practicable because the statistics (with the standard segment density of 10) are not sufficient. It so appears that the results obtained so far are valid for well entangled polymers only, which is indeed the main aim of the present contribution. Finally, Fig. 5 shows the kind of microphase structures obtained, depending on block ratio. Reasonable reproductions of the sphere and cylinder phases can be recognized. However, in the present work the size of the simulation box was fixed independently of the equilibrium microphase structure to be eventually reached, and was not optimized.
80 SCF N = 40n0 N = 20n0
60
χN
Hence the ordered phases, as obtained in classical density field calculations, were not exactly reproduced. 5. Discussion
S
C
L
40 G 20
0
Fig. 5. Snapshots after 100 000 time units for diblock copolymers with the same molecular mass as that in Fig. 1 (N = 40n0), but with different block ratios and interaction parameters: (a) / = 0.1 and vn0 = 2.0; (b) / = 0.3 and vn0 = 2.0.
DISORDERED 0.0
0.1
0.2 0.3 block ratio
0.4
0.5
Fig. 4. The phase diagram for two diblock copolymers with N = 40n0 and N = 20n0. Filled symbols indicate microphase separation and open ones correspond to the disordered state. SCF results [23] are also presented for comparison (lines) with the indication of the ordered phase structure: (S) sphere, (C) cylinder, (L) lamella, and (G) gyroid.
In this section, we discuss how the proposed model compares with, and possibly complement, the conventional methods employed for block copolymers. It is emphasized that our method is the single molecularly-based model able to describe with acceptable computational cost the global motion of well-entangled block copolymers. Indeed computations were performed with a Linux workstation equipped with a 3.0 GHz Pentium 4 processor with 1 GB main memory. With that computer, and for the 163 simulation box filled with the typical segment density of 10 (i.e., for a total number of segments of 10 · 163 = 4.1 · 104), ca. 16 h are required for a run of 5000 time units, roughly corresponding to the longest relaxation time of chains with 40 entanglements. It is also noteworthy that the calculation time is proportional to
Y. Masubuchi et al. / Journal of Non-Crystalline Solids 352 (2006) 5001–5007
the system size since the model does not account for twobody non-bond interactions, which are typically the most expensive in conventional molecular dynamics as their number increases with the square of the system size. However, as mentioned before, our computations have shown some difficulties with marginally entangled chains because of the resolution of the chemical potential field. Also in view of the fact that our model certainly does not apply to unentangled chains, multi-scale simulations are required if one needs to explore the whole range of molecular weights. Of course, our computations of the dynamical evolution are not fast if compared with the density field calculations used to obtain equilibrium structures. For instance, a recently reported simulation [24] takes 3 h (with a similar computer to that used in the present work) to obtain the equilibrium structure of a mixture of block-copolymers and linear polymers in a box of 643 lattice points. On the other hand, since most density-field-based dynamic simulations do not account for entanglements, a suitable combination of different methods complementing each other seems promising. Indeed, the equilibrium structure quickly obtained with the density field theories certainly remains useful because the present method cannot easily reproduce the exact ordered structures of the microphases. 6. Conclusions A simulation method for entangled block copolymers was here proposed as an extension of the primitive chain network model [10,15]. Predictions for the phase diagram are in qualitative agreement with the conventional methods, but phase separation kinetics is also described. Hence, the present method may occupy an important niche among the existing ones, falling somewhere in between molecular dynamics [6], density field calculations [2], and singlechain sliplink simulations [9,11]. Detailed assessment of the method by quantitative comparison with experiments on rheology, phase separation kinetics, molecular confor-
5007
mation and molecular motion, etc. is in progress, and will be presented elsewhere. Acknowledgement This work was supported in part by a Grant in Aid for Scientific Research [KAKENHI] 2005 No. 17750200 from the Ministry of Education, Culture, Sports, Science and Technology of Japan. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24]
G. Srinivas, D.E. Discher, M.L. Klein, Nature Mater. 3 (2005) 638. G.H. Fredrickson, V. Ganesan, F. Drolet, Macromolecules 35 (2002) 16. K. Binder, Comp. Phys. Commun. 147 (2002) 22. T. Dotera, A. Hatano, J. Chem. Phys. 105 (1996) 8413. A. Terzis, D.N. Theodorou, A. Stroeks, Macromolecules 35 (2002) 508. M. Murat, G.S. Grest, K. Kremer, Macromolecules 32 (1999) 595. P. Ahlrichs, B. Dunweg, J. Chem. Phys. 111 (1999) 8225. T. Aoyagi, T. Honda, M. Doi, J. Chem. Phys. 117 (2002) 8153. C.C. Hua, J. Schieber, J. Chem. Phys. 109 (2001) 10018. Y. Masubuchi, J.-I. Takimoto, K. Koyama, G. Ianniruberto, G. Marrucci, F. Greco, J. Chem. Phys. 115 (2001) 4387. J.-I. Takimoto, M. Doi, Philos. Trans. Royal Soc. London A 361 (2003) 641. J.D. Schieber, J. Neergaard, S. Gupta, J. Rheol. 47 (2003) 213. A.E. Likhtman, Macromolecules 38 (2005) 6128. M. Doi, S.F. Edwards, The Theory of Polymer Dynamics, Clarendon, Oxford, 1986. Y. Masubuchi, G. Ianniruberto, F. Greco, G. Marrucci, Model. Simul. Mater. Sci. Eng. 12 (2004) S91. T. Shima, H. Kuni, Y. Okabe, M. Doi, X.F. Yuan, T. Kawakatsu, Macromolecules 36 (2003) 9199. J. des Cloizeaux, Europhys. Lett. 5 (1988) 437. G. Marrucci, J. Non-Newtonian Fluid Mech. 62 (1996) 279. J.T. Padding, W.J. Briels, J. Chem. Phys. 117 (2002) 925. P.J. Flory, Principles of Polymer Chemistry, Cornell University Press, Ithaca, 1953. Y. Masubuchi, G. Ianniruberto, F. Greco, G. Marrucci, KoreaAustralia Rheology J. 18 (2006) 99. T. Yaoita, T. Isaki, Y. Masubuchi, H. Watanabe, G. Ianniruberto, F. Greco, G. Marrucci, J. Chem. Phys. 121 (2004) 12650. M.W. Matsen, F.S. Bates, Macromolecules 29 (1996) 1091. T. Uneyama, M. Doi, Macromolecules 38 (2005) 196.