J. Phys. Chem. B 2007, 111, 13937-13942
13937
Mesoscopic Simulation on Phase Behavior of Pluronic P123 Aqueous Solution Yurong Zhao, Xiao Chen,* Chunjie Yang, and Guodong Zhang Key Laboratory of Colloid and Interface Chemistry, (Shandong UniVersity), Ministry of Education, Jinan, Shandong, 250100, Peoples Republic of China ReceiVed: May 6, 2007; In Final Form: October 11, 2007
As a new dynamic density functional method, the mesoscopic dynamics (MesoDyn) is used to simulate the microphase separation of the binary mixture of tri-block copolymer P123 ((EO)20(PO)70(EO)20) and water. With a simple copolymer model, various aggregate structures of P123 in water including the micelle, hexagonal, and lamellar phases are produced, which can partly reproduce most experimental phase regions. The simulated phase ranges are more or less different from those established from experiment, especially at high polymer concentrations. This can be attributed to different phase mapping situations, that is, the constant shear used in simulation versus the varied external forces in experiment. The calculated trend for micelle size change in diluted region is co-incident well with previous observations in other Pluronic systems. The poly(propylene oxide) (PO) block amount is found to influence not only the aggregate morphology but also its formation rate in solution. Obtained results indicate that the mesoscopic simulation is a valuable tool to supplement the experimental study on aggregates formation.
1. Introduction Pluronic P123 ((EO)20(PO)70(EO)20) [block PEO-PPO-PEO copolymer, PEO ) poly(ethylene oxide), PPO ) poly(propylene oxide)] is an important non-ionic surfactant that has been widely used in detergency, foam formation, dispersion stabilization, lubrication, and drug delivery.1-5 Depending on the temperature and concentration, it is capable of building diverse structures in water as can be seen from a typical phase diagram of P123water system determined by Wanka et al.6 With P123 concentration increasing, different aggregates such as micelle, cubic, hexagonal, and lamellar phases are formed in this binary system. A variety of experimental techniques such as small-angle neutron scattering, small-angle X-ray scattering, static or dynamic light scattering, rheology, polarizing optical microscopy, differential scanning calorimetry, 1H NMR spectrum, and surface tension measurements have been used to investigate the structural and dynamical properties of P123-water system.7 Though obtained results have obviously improved the understanding for the phase behavior of this binary system, its intrinsic microphase separation information is still ambiguous, and the theoretical simulation research is also rarely reported. Recently, with the fast development of computers, simulation has become a powerful complementary tool for the study of surfactant systems at a molecular level, like molecular dynamics (MD) simulations.8-12 However, building a surfactant mesophase model with classical MD methods at atomic resolution is currently not possible because of the time and length scale limits at which these complex phase structures can be formed. By taking into account all of these difficulties, an alternative and effective mesoscopic dynamics technique, MesoDyn, is more adopted to simulate the phase behavior of aqueous Pluronic solutions. On the basis of the dynamic mean-field density function theory and the characteristic merit to run the time integration of functional Langevin equations, MesoDyn simula* To whom correspondence should be addressed. E-mail: xchen@ sdu.edu.cn. Fax: +86-531-88564464. Tel: +86-531-88365420.
tion has led to significant advances in the investigation for the microphase separation of block copolymers,13-17 including simulations on the time evolution of micelle or mesoscale structure formation, phase behavior of block copolymers, drug delivery, and the solution behavior of specific compounds.14-17 In contrast to previous simulation approaches aimed at classifying morphologies by means of equilibrium theories, the MesoDyn method recognizes the fact that these patterns are irregular in nature and, hence, can only be characterized via the dynamic properties of the systems. This approach is much more reliable from an industrial perspective, since typical processing times are orders of magnitude shorter than the thermodynamic relaxation time, and such nonperfect states thus contribute substantially to the behavior of the final material.18 Another important advantage of this method is that no prior assumption is needed on the phases, and the kinetics of phase formation can thus be studied, which is very difficult to observe in experiment. The aim of our present MesoDyn study is to simulate the phase behavior of a binary system comprising P123 and water, which has been done from the homogeneous density distribution by an instantaneous quench. Much interest is focused on the existing shear at high concentrations. Influences of polymer concentration or PO block content on the aggregate morphology and formation rate are also discussed. It is believed that such simulation will provide a complementary way for better understanding the phase behavior of this binary system and the related microphase separation information. 2. Simulation Method and Details 2.1. Theory in MesoDyn. The basic idea in the MesoDyn method is the density function theory, which is based on the concept that the free energy F of an inhomogeneous liquid is a function of the local density function F. From the free energy, all thermodynamic functions can be derived. The model used in the MesoDyn project consists of beads of various types I, J, ... with interactions described by harmonic
10.1021/jp073451f CCC: $37.00 © 2007 American Chemical Society Published on Web 11/29/2007
13938 J. Phys. Chem. B, Vol. 111, No. 50, 2007
Zhao et al.
Figure 2. Density slices respectively for PEO, PPO, and water in micelles at low concentrations.
relation between densities, distributions, and external potential, can be written as: Figure 1. Simulated phase evolution of Pluronic P123 with its concentration (in volume fraction) in aqueous solution at 298 K. M, H, and L denote regions respectively for micelle, hexagonal, and lamellar phases. Shear is applied at the concentration >30%.
oscillator potentials for intramolecular interactions (Gaussian chain) and a mean field potential for all other interactions.19 Each bead is of a certain component type representing covalently bonded groups of atoms such as those given by one or a few structural units of a polymer chain. The dynamics of the system is described by a set of functional Langevin equations, which, in simple terms, are diffusion equations in the component densities by taking account of the noise in the system. On a coarse-grained time scale, F0(r) is defined as a collective concentration field for the type I beads at an instant of time and serves as a reference level. There will be a certain distribution of bead positions, defined as Ψ(R11,...,RnN) where Rγs is the position of bead s from chain γ. Given the distribution Ψ, the collective concentration of the bead s from all chains can be defined by the average of a microscopic density operator: n
F1[Ψ](r) ≡
N
δK1STrΨδ(r - Rγs) ∑ ∑ r)1 s)1
(1)
where δK1S is the Kronecker function with value 1 when bead s is of type I and 0 otherwise. It is assumed that in the liquid with slow relaxation the interactions do not depend on the momenta. A set of distribution functions Ψ is defined with the constraint F01(r) ) F1[Ψ](r). All distributions Ψ lead to the same density F01(r) to form an equivalent class Ω of distribution functions:
Ω ) {Ψ(R11, ..., RnN)|F1[Ψ](r) ) F01(r)}
β-1 ln n! -
(3)
The first term is the average value of the Hamiltonian for internal Gaussian chain interactions.20,21 The second and third terms represent respectively the Gibbs entropy of the distribution -kBTΨ ln Ψ and the mean-field nonideal contribution. Ψ is independent of the system history and is fully characterized by the constraint that represents the density distribution and minimizes the free-energy function. This constraint on the density fields is realized by means of an external potential UI. The constraint minimization of the free energy function leads to an optimal distribution, which in turn, by the one-to-one
∑1 ∫U1(r)F1(r) dr + βFnid[F]
(4)
Now the nonideal free energy function Fnid is introduced,
∫∫
1 ∈AA(|r - r′|)FA(r)FA(r′) + 2 ∈AB(|r - r′|)FA(r)FB(r′) + ∈BA(|r - r′|)FB(r)FA(r′) + ∈BB(|r - r′|)FB(r)FB(r′) dr dr′ (5)
Fnid[F] )
where ∈IJ(|r - r′|) is a mean-field energetic interaction between beads of type I at r and type J at r′. The mean-field intrinsic chemical potentials can easily be derived by functional differentiation of the free energy µI(r) ) δF/δF1(r). At equilibrium, µI(r) ) constant, which results in the familiar self-consistent field equations for the mean-field Gaussian chain model. In general, these equations have many solutions, one of which will be a state of lowest free energy, while most other states will be metastable. On the basis of these equations, the generalized time-dependent Ginzburg-Landau theory can be set up. The following functional Langevin equations for the diffusive dynamics of the density fields are also introduced:
∂FA ) MυB∇FAFB∇[µA - µB] + η ∂t
(6)
∂FB ) MυB∇FAFB∇[µB - µA] + η ∂t
(7)
The distribution of the Gaussian noise η satisfies the fluctuation-dissipation theorem:
(2)
On the basis of this set of distribution functions, an intrinsic free-energy function F[Ψ] can be defined:
F[Ψ] ) Tr(ΨHid + β-1Ψ ln Ψ) + Fnid[F0]
βF[F] ) n ln Φ +
〈η(r,t)〉 ) 0 〈η(r,t)η(r′,t′)〉 ) -
(8)
2MυB δ(t - t′)∇r × δ(r - r′)FAFB∇r′ (9) β
where M is a bead mobility parameter. The kinetic coefficient MνFΑFΒ models a local exchange mechanism. The Langevin equations are constructed for an incompressible system with dynamic constraint:
(FA(r,t) + FB(r,t)) )
1 υB
(10)
where υB is the average bead volume. 2.2. Simulation Parameters. In MesoDyn simulations, two sets of parameters must be defined to specify the chemical nature of the system: one is the chain topology in terms of repeat
MesoDyn Simulation of P123 Aqueous Solution
J. Phys. Chem. B, Vol. 111, No. 50, 2007 13939
Figure 3. Time evolution of dimensionless order parameters for P123-water system at different concentrations: (a) 10% and (b) 20%. Simple kinetics models formed by 10% P123-water system at different time steps (i)100, (ii)1000, and (iii) 14 000 are given below the curves; these structures correspond to the stages I, II, and III, respectively.
Figure 4. Different morphologies at variable shear rates (s-1) for 45% P123-water system at the same condition: (a) 0, (b) 1 × 102, (c) 1 × 103, (d) 1 × 104, (e) 1 × 105, and (f) 1 × 106.
segments (or beads), and the other is the interaction energy between different components. For the first set, MesoDyn uses a Gaussian chain “springs and beads” description, where each bead in the Gaussian chain is a statistical unit, representing a number of “real” monomers and different types of beads correspond to diverse components. Because springs reflect the stretching behavior of a chain fragment and the chain topology depends on the coarsening degree of the original system, the choice of the Gaussian chain becomes much more important. van Vlimmeren et al.22 established a simple relationship between the atomic chains and Gaussian chains for Pluronics, which gave the numerical values corresponding to the number of the units/ bead ratio of 4.3 for two PEO (poly(ethylene oxide)) blocks and 3.3 for a PPO (poly(propylene oxide)) block. The numbers of PEO/PPO units in P123 or F127 are 20/70 or 100/70 respectively. Thus, we can get their coarsened chain topology numbers of 4.7/21.2 for P123 and 23.3/21.2 for F127. However,
because both numbers should be integers, we choose the Gaussian chains respectively as A5B21A5 for P123 and A23B21A23 for F127, where the solvophilic A and solvophobic B blocks represent PEO and PPO segments. Moreover, one bead W represents water in MesoDyn simulations. Such choice of Gaussian chains in our study is confirmed reasonable if comparing the simulated morphology to the experimental results.6,23 The interaction energy ∈IJ (in eq 5) between different kinds of segments in MesoDyn represents the pairwise interactions of beads, which is similar to that defined in the Flory-Huggins model.22 It can be derived either from atomic simulation,24,25 empirical methods,26 or from experimental data such as vapor pressure.22 The simplest method is based on regular solution theory and relates the Flory-Huggins parameter χ to the component solubility parameter δ. These parameters can be considered as the nonideal interactions that are included via a mean-field approach, and the strength of repulsion interaction between different components is characterized by >0 in units of kiloJoules per mole. The solvent-polymer interaction parameter is estimated as:
χsp ) (δs - δp)2Vs/RT
(11)
where δs and δp are respectively Scatchard-Hildebrand solubility parameters of the solvent and the polymer.27 Vs is the molar volume of the solvent. In our study, the interaction parameters between H2O, PPO, and PEO are not calculated according to eq 11. They are directly selected as the same values as those adopted in previous successful simulations for the Pluroric solutions.13-17,22,24 The Flory-Huggins parameters are chosen as χAW ) 1.4, χBW ) 1.7 for solvent-polymer interaction and χAB ) 3.0 for polymerpolymer interaction. For all simulations, the dimensionless
13940 J. Phys. Chem. B, Vol. 111, No. 50, 2007
Zhao et al.
Figure 5. PPO isosurface representations of P123-water system after 20 000 simulation steps in the absence of shear: (A) 7%, (B) 10%, and (C) 20%.
parameters in MesoDyn program have been chosen as the time step ∆τ ) 50 ns (dimensionless time step ∆τ ) 0.5), the noise scaling parameter Ω ) 100, and the compressibility parameter κ′H ) 10; the grid parameter d is automatically set to 1.15430 to ensure isotropy of all grid-restricted operators,13 and the total simulation time is 5000 steps with shear and 20 000 steps without shear. The simulation temperature is set to 298 K, and the bead diffusion coefficient D ) β-1M of the three types is set to 1.0 × 10-7cm2 s-1. All simulations are processed in a cubic grid with 32 × 32 × 32 cells of mesh size h and carried out with MesoDyn module in a commercial software package Cerius2, version 4.6, from Accelrys, Inc.28 3. Results and Discussion 3.1. Phase Behavior of Aqueous Pluronic P123 Solution. It is well-known that aqueous Pluronic P123 solution is less viscous at low concentrations, and it is comparatively easy to make such systems homogeneous without any external force or with only a little one. At higher concentrations, for example, above 30% (volume fraction, thereinafter), however, the solution will become more and more viscous and an external force must be applied to make it homogeneous. With these in mind, a shear (γ ) 1 × 106 s-1) is added at higher concentrations (>30%) in our MesoDyn approach, which is equivalent to the external force used in experiment. Figure 1 shows the simulated mesoscale structures as a function of block-copolymer concentration. It is clearly seen that different aggregates such as micelle, hexagonal, or lamellar phases are gradually formed in aqueous solutions. From the density slices of PEO, PPO, and water in micelles (see Figure 2), it is easy to find that the hydrophobic PO blocks are gathered into a spherical core, while the hydrophilic EO blocks are solvated by water, exhibiting most characteristics of spherical micelles. But, some experiment results have certified the existence of anisometric micelles with an axial ratio slightly larger than 1 in the P123-water system,6 which may be due to the larger aggregation number of P123 (about 3 times higher than those of other Pluronic polymers). Such anisometric micelles can be thought of as one type of irregular spherical ones, pseudo-spherical micelles, which may have some different properties from the conventional spherical micelles. In order to investigate the simulated micelles deeply, we calculated the aggregation number of these micelles. In a 32 nm × 32 nm × 32 nm cubic grid, the number of the micelle is averaged, about 20, for 3-29 vol % systems. Then, the average core radius, Rc, of PO monomer is calculated, which is about Rc ) 73 Å. To a good approximation, the aggregation number (N) scales linearly with the third power of the core radius, Rc:29
4 NnVPO ) πRc3 3
(12)
where n ) 70 is the PO degree of polymerization, and the propylene oxide volume is VPO ) 95.4 Å3. This leads to an
aggregation number of N ≈ 245 for the micelles. It is very close to the value determined by static light scattering measurements at 300 K.6 This calculation of the aggregation number helps us better understand these particular micelles. After the micelle, the hexagonal or lamellar phases are subsequently formed respectively by parallelly arranged long rod-like micelles or the stacked bilayers. It can be seen that these three different phases (micelle, hexagonal, and lamellar) occur respectively at various concentration regions, like 3∼29% for pseudo-spherical micelle and 32∼55% for hexagonal phase. Though the calculation denotes 65% as the beginning polymer concentration for the lamellar phase, a regular one can still not be simulated until 80%. Compared with the corresponding experimental phase diagram,6 our simulation can partly reproduce most phase regions, except for the absence of cubic phase. The simulated phase regions are more or less different from those established by experiment, especially at high concentrations, which can be attributed to different phase mapping situations, that is, the constant shear used in simulation versus the varied external forces in experiment. In order to prove the effect of the shear rate, we calculated most of our systems with variable shear rates. Simulation results show that different phases can be obtained at various shear rates for the same concentration. Now, we set a 45% P123-water system as an example to illuminate the influence of variable shear rates to the phase behavior of the Pluronics system. To carry this out, we add different shear rates for this system, from 0 s-1 to 1 × 106 s-1. Figure 4 shows different morphologies at variable shear rates. It is easy to find that, when the shear rate is 0, bicontinuous phases can be formed; upon increasing the shear rate, the morphology changes, but the essential structure keeps. When the shear rate arrives at 1 × 105 s-1, the regular hexagonal phase is formed. To keep on increasing the shear rate, the essential hexagonal phase stays. From this, why experimental and simulation results have differences can be understood to a certain extent. Besides, chain model and interaction parameters may also be very important causes of discrepancies between simulation and experiment. What’s more, hydrodynamic interaction may also be a possible reason, but as for our symmetric polymer, the effect may be rather little.30 As a whole, though the theoretical method may have some limitations, it can still be an effective supplement for experiments. We will do much work to improve this in the future. 3.2. Phase Evolution and the Concentration Effect. The order parameter P, representing the characteristics of phase separation, is defined as follows:
∫V ∑ [θ21(r) - (θ01)2] dr P≡
1
V
(13)
where θ is the polymer volume fraction and V is the volume of cells. P is the mean-squared deviation from homogeneity in the system, which captures both the effects of phase separation and compressibility. In order to show the process of phase separation at different concentrations more distinctly, it is helpful to analyze the order parameter changes during the formation of the aggregates. The time evolutions of order parameters for 10% and 20% P123 aqueous solutions are shown in Figure 3, which exhibits three stages (I, II, and III in Figure 3a) in the phase evolvement. Pre-micelles are formed in the first stage (I), where the order parameter and the aggregate morphology change little, corre-
MesoDyn Simulation of P123 Aqueous Solution
J. Phys. Chem. B, Vol. 111, No. 50, 2007 13941
Figure 6. Time evolution of the dimensionless order parameters of two systems under the same conditions but with different PO content: (a) 10% F127-H2O and (b) 10% P123-H2O.
sponding to the nucleation of several polymer aggregates. Then, the order parameter increases rapidly, indicating the formation of micelles at the second stage (II), showing their approximate morphology. The period of this stage is about tens of microseconds. The final stage (III) is time-consuming, at which the system changes slowly to overcome the defects formed in the previous stage. Simple kinetics models formed by a 10% P123water system at these three different time steps, 100, 1000, and 14 000, are given below the curves; these structures correspond to the stages I, II, and III, respectively. Now, these three stages will be clear at a glance. For systems with the hexagonal or lamellar phase formation, the order parameters go through in a similar way. By comparing Figure 3a,b, it is easy to find that the higher polymer concentration will make the phase separation much faster because of the stronger thermodynamic driving forces.17 This indicates that the polymer concentration affects not only the presence of different phases but also their formation rates, which is difficult to observe directly by experiment. 3.3. Effect of P123 Concentration on Micelle Size. For P123, the concentration effect on the micelle core size is simulated in the region M shown in Figure 1. Such core size is found to increase with the polymer concentration (see Figure 5), which is considered to be caused by a concomitant decrease of the water content and therefore a less strong poor solvent effect on the PPO block core. The core is thus not going to collapse as much as those at lower polymer concentrations,24 which will induce bigger cores and result in larger micelles. This result is shown to compare well with previous observations on F127 [(EO)100 (PO)70(EO)100] aqueous solution by both theoretical24 and experimental31 techniques. On the other hand, Wanka et al.6 deduced on the basis of SANS data that the size of micelles is independent of concentration. In SANS experiments, a hard sphere is used to fit for the form factor; that is, a sharp interface between PEO and PPO is supposed. The discrepancy from our simulation arises from the fact that the fitting of SANS data ignores the hydrated layer of PEO at the interface between the hydrophobic core (PPO chains) and hydrophilic corona (PEO chains).32 Therefore, Mesodyn could provide results corresponding to relatively more actual situations. 3.4. Effect of PEO and PPO Relative Block Size. In order to investigate the effect of relative sizes of the hydrophilic EO
block and hydrophobic PO block on the phase behavior (the formation of ordered aggregates), a new system (10% F127 aqueous solution) under the same conditions but with different block sizes is simulated. From the dimensionless order parameters of two systems, that is, 10% F127 and 10% P123 aqueous solutions (see Figure 6a,b respectively), it is easy to find that, with the less hydrophobic PO component, F127 is more difficult to form ordered aggregates and makes the order parameter get equilibrated in a longer time. This is because the decrease of PO block content is accompanied correspondingly by a decrease of the hydrophobic counterpart and hence the enhanced difficulty to form the hydrophobic core and the phase separation. It is obvious that the order parameters from simulation can supply a much easier and more intuitionistic way to understand the effect of PEO and PPO relative block size on the speed of aggregate formation. 4. Conclusions The phase behavior of a Pluronic P123-water system in the mesoscopic region is successfully simulated using the MesoDyn equivalent chain method. Obtained results show that, with the increase of concentration of P123, different aggregates such as micelle, hexagonal, and lamellar phases are formed in aqueous solutions. The influences of PO block amount or the relative amounts between PO and EO blocks on aggregate morphology and formation rate can be reasonably simulated. The simulated results are more or less different from those established by experiment, especially at high concentrations, reminding us that the actual systems are difficult to reach the real thermodynamic equilibrium and shears should be considered in simulation. It can be concluded that the mesoscopic simulation method is a valuable tool for the description of mesoscale morphology formation and gives an insight into the process of aggregate formation. Acknowledgment. We are thankful for the financial supports from the National Natural Science Foundation (20373035, 20573066) and Science Fund of Shandong Province (Y2005B18) in China. References and Notes (1) Tuzar, Z.; Kratochvil, P. AdV. Colloid Interface Sci. 1976, 6, 201. (2) Zhou, Z.; Chu, B. Macromolecules 1994, 27, 2025.
13942 J. Phys. Chem. B, Vol. 111, No. 50, 2007 (3) Chu, B.; Zhou, Z.; Wu, G.; Farrell, H. M. J. Colloid Interface Sci. 1995, 170, 102. (4) Brown, W.; Schillen, K.; Almgren, M.; Hvidt, S.; Bahadur, P. J. Phys. Chem. 1991, 95, 1850. (5) Wanka, G.; Hoffmann, H.; Ulbricht, W. Colloid Polym. Sci. 1990, 268, 101. (6) Wanka, G.; Hoffmann, H.; Ulbricht, W. Macromolecules 1994, 27, 4145. (7) Ian W. Hamley, Block Copolymers in Solution: Fundamentals and Applications; John Wiley & Sons: Canada, 2005. (8) Alaimo, M. H.; Kumosinski, T. F. Langmuir 1997, 13, 2007. (9) Allen, R.; Bandyopadhyay, S.; Klein, M. L. Langmuir 2000, 16, 10547. (10) Bandyopadhyay, S.; Tarek, M.; Lynch, M. L.; Klein, M. L. Langmuir 2000, 16, 942. (11) Salaniwal, S.; Cui, S. T.; Cochran, H. D.; Cummings, P. T. Langmuir 2001, 17, 1784. (12) Abel, S.; Sterpone, F.; Bandyopadhyay, S.; Marchi, M. J. Phys. Chem. B 2004, 108, 19458. (13) Fraaije, J. G. E. M.; van Vlimmeren, B. A. C.; Maurits, N. M.; Postma, M.; Evers, O. A.; Hoffmann, C.; Altevogt, P.; Goldbeck-Wood, G. J. Chem. Phys. 1997, 106, 4260. (14) Fraaije, J. G. E. M.; Sevink, G. J. A. Macromolecules 2003, 36, 7891. (15) Lyakhova, K. S.; Zvelindovsky, A. V.; Sevink, G. J. A.; Fraaije, J. G. E. M. J. Chem. Phys. 2003, 118, 8456. (16) Altevogt, P.; Evers, O. A.; Fraaije, J. G. E. M.; Maurits, N. M.; van Vlimmeren, B. A. C. J. Mol. Struct. (THEOCHEM) 1999, 463, 139.
Zhao et al. (17) Zhang, X. Q.; Yuan, S. L.; Wu, J. Macromolecules 2006, 39, 6631. (18) Li, Y. M.; Xu, G. Y.; Chen, A. M.; Yuan, S. L.; Cao, X. R. J. Phys. Chem. B 2005, 109, 22290. (19) Doi, M.; Edwards, S. F. The Theory of Polymer Dynamics; Oxford Science Publications: Oxford, United Kingdom, 1986. (20) Barber, M. N.; Ninham, B. W. Random and Restricted Walks; Gordon and Breach: New York, 1970. (21) Abe, A.; Mark, J. E. J. Am. Chem. Soc. 1976, 98, 6468. (22) van Vlimmeren, B. A. C.; Maurits, N. M.; Zvelindovsky, A. V.; Sevink, G. J. A.; Fraaije, J. G. E. M. Macromolecules 1999, 32, 646. (23) Holmqvist, P.; Alexandridis, P.; Lindman, B. J. Phys. Chem. B 1998, 102, 1149. (24) Lam, Y. M.; Goldbeck-Wood, G. Polymer 2003, 44, 3593. (25) Zhang, M.; Choi, P.; Sundararaj, U. Polymer 2003, 44, 1979. (26) Honeycutt, J. D. Comput. Theor. Polym. Sci. 1998, 8, 1. (27) Barton, A. F. M. Handbook of Solubility Parameters and Other Cohesion Parameters; CRC Press: Boca Raton, FL, 1983. (28) One Molecular Simulation Software, Inc., see https://www.accelrys.com. (29) Mortensen, K.; Pedersen, J. S. Macromolecules 1993, 26, 805. (30) Groot, R. G.; Maden, T. J.; Tildesley, D. J. J. Chem. Phys. 1999, 110, 9739. (31) Lam, Y. M.; Grigorieff, N.; Goldbeck-Wood, G. Phys. Chem. Chem. Phys. 1999, 1, 3331. (32) Mortensen, K. J Phys: Condens. Matter 1996, 8, 103.
