Diffusion of Nanoparticles in Semidilute Polymer Solutions:A Multiparticle Collision Dynamics Study

Shu-xian Li,Hui-jun Jiang,Zhong-huai HouDepartment of Chemical Physics and Hefei National Laboratory for Physical Sciences at the Microscales, University of Science and Technology of China,Hefei 230026,China

Diffusion of Nanoparticles in Semidilute Polymer Solutions:A Multiparticle Collision Dynamics Study

Shu-xian Li,Hui-jun Jiang,Zhong-huai Hou∗
Department of Chemical Physics and Hefei National Laboratory for Physical Sciences at the Microscales, University of Science and Technology of China,Hefei 230026,China

The diffusion of nanoparticles immersed in semidilute polymer solutions is investigated by a hybrid mesoscopic multiparticle collision dynamics method.Effects of polymer concentration and hydrodynamic interactions among polymer monomers are focused.Extensive simulations show that the dependence of (diffusi)on coefficient D on the polymer concentration c agrees with Phillies equation D-exp−αcδwith a scaling exponent δ≈0.97 which coincides with the experimental one in literature.For increasing nanoparticle size,the scaling prefactor α increases monotonically while the scaling exponent always keeps fixed.Moreover,we also study the diffusion of nanoparticle without hydrodynamic interactions and find that mobility of the nanoparticle slows down,and the scaling exponent is obviously different from the one in experiments,implying that hydrodynamic interactions play a crucial role in the diffusion of a nanoparticle in semidilute polymer solutions.

Nanoparticle,Polymer solution,Multiparticle collision dynamics

I.INTRODUCTION

Understanding the transport properties of nanoparticles(NPs)in polymer liquids,including polymer solutions and polymer melts,is of great importance in many interdisciplinary fields.For instance,the diffusive behavior of NPs can provide important information about the local mechanical and viscoelastic properties of polymer liquids[1-3].Diffusion dynamics of a probe in dense polymer solutions can help to understand effects of crowding on intracellular diffusion processes of proteins or other macromolecules[4-7].Adding NPs to polymer liquids can result in novel electrical and photonic properties of nanocomposite systems where the mobility of NPs can play a crucial role[8-12],to list just a few.

For many of these reasons,diffusion of NPs in polymer solutions has received a lot of attention both experimentally[13-21]and theoretically[22-26].It was found that the well-known Stokes-Einstein(SE)relationship between the long-time diffusion coefficient D and the particle size might be violated in polymer solutions and be strongly dependent on the interactions between the probe particle and the polymer molecules[13,19].In particular,much attention has been paid to how D depends on the concentration c of polymer solutions by using experimental techniques such as fluorescence correlation spectroscopy measurements.In the semi-dilute transition regime,the dependence of D on c was shown to follow a stretched exponen(tial fu)nction given by the Phillies equation[27],D-exp−αcδ,the exponent δ is usually less than one and α is dependent on the probe size[28].Such an observation was also confirmed by a recent experiment[16],the diffusion of gold NPs of radii 2.5−10 nm in semidilute poly(ethylene glycol)(PEG) water solution was measured.Theoretically,NP's diffusion in polymer liquids has been studied by using scaling theory[22,26,29],“walking confined diffusion”model [23,24],mode coupling theory(MCT)[25,30-33],and so on.In the scaling theory developed by Cai et al. [22],the effect of chain relaxation on the mobility of nonsticky NPs in polymer liquids was considered,and a power law dependence of the diffusion coefficient on polymer concentration was derived.The“walking confined diffusion'‘'model[23,24]considered a depletion layer and it could reproduce the nonlinear dependence of mean square displacement on time.In recent years, an important theoretical framework based on MCT has been proposed to study the long-time diffusion coefficient in complex polymer liquids,from a microscopic level explicitly accounting for the solute-solvent interactions.It starts from the calculation of the friction kernel ζ(t),defined as the time correlation function of the random force exerted on NP,which may result from short time collisions among NP and surrounding particles,the coupling of NP dynamics with the solvent structural relaxation modes,or coupling with a long-time transverse hydrodynamic currents.With the calculated ζ(t), the diffusion coefficient D can be obtained directly viafluctuation-dissipation theorem.Using MCT,Schweizer and coworkers[25,31]had studied the diffusion of NPs in polymer melts theoretically,being both unentangled and entangled,which helps to understand many important experimental observations and also provides deep physical insights.Egorov[30]performed a slightly different MCT study on the anomalous diffusion behavior of NP in polymer melts and solutions,but only unentangled regime was considered and the effect of solvent was implicitly accounted for.Most recently,Dong et al.extended the MCT to calculate D of NPs in PEG solutions by introducing a particular dynamic scattering function suitable for polymer solutions,finding very good quantitative agreements with experimental data [33].In parallel to the theoretical approaches,there have been some coarse-grained molecular dynamics simulations(CGMD)as well[11,34-37].Liu and his coauthors found that the gyration radius of a polymer chain is a key factor to determine the validity of the SE relationship in describing NPs diffusion in polymer melts [35].While a recent CGMD study performed by Kalathi et al.showed that it is the entanglement mesh size that should be such a key factor and an NP's diffusivity has two very different classes of behavior depending on its size[37].

∗Author to whom correspondence should be addressed.E-mail: hzhlj@ustc.edu.cn

We note that diffusion of NPs in polymer solutions is a rather complicated problem,wherein not only the interactions between NPs and polymer molecules must be considered,but also the solvent effects may play important roles.Nevertheless,most of the theoretical and simulative studies so far have not related to solvent effects explicitly.MCT studies mentioned above, only considered either polymer melts or the solvent implicitly.And to the best of our knowledge,studies on the diffusion behavior of NPs by using CGMD method mainly concerned about polymer melts rather than solutions.On the other hand,it has been shown that solvent effects,particularly the long-range hydrodynamic interaction(HI),might drastically influence dynamics of polymer chain in a solvent[38-42].For example, Kikuchi et al.reported that HIs accelerate the coilglobule transition(CGT)of a flexible polymer and also affect the morphology in the shrinking process[38,39]. Chang and Yethiraj also proposed that HIs tend to prevent a collapsing polymer from being trapped at local energy minima[40].Recently,Tanaka speculated that HIs not only accelerate the CGT,but also may retard it via a squeezing flow effect[41].Moreover,Kamatu and coworkers found that HIs accelerate collapsing for a quench from above Flory temperature(i.e.θ point),whereas they decelerate collapsing for a quench from below θ point and they believed that the roles of HIs in the chain collapsing transition crucially depend upon the initial enhancement of anisotropy of a polymer configuration[42].These studies clearly showed HIs can play subtle roles in both polymer collapse and the folding pathway.Thus,it is quite natural for one to seek for some mesoscopic simulation approaches to study dynamics of such complex systems.Recently,a new approach termed as multiparticle collision dynamics(MPCD)was proposed to address this issue.The main idea of MPCD is to replace the solvent molecules by coarse-grained particles,which performs streaming and collision steps to make sure the long range HI effects are maintained[43,44].This mesoscopic method, sometimes combined with MD,has been used successfully in many soft matter systems[45,46],including those involving polymers[44,47,48].

In the present work,we investigate the diffusion of an NP in a semi-dilute polymer solution by using the MPCD method combined with MD.The main motivation of the present study is to take into account the effects of solvent in an explicit and direct way,including the particular role of HI.For simplicity,we only consider the semi-dilute regime and mainly investigate the dependence of diffusion coefficient D on the polymer concentration.We find that our method can well reproduce the scaling relation between the gyration radius and the concentration.The diffusion coefficient, as a function of polymer concentration c,is shown to follow the Phillies equation very well with reasonable exponents α and δ in agreements with experimental observations.Importantly,the method facilitates us to study the very role of HI by switching on or off the HI during the simulation process.It is demonstrated that D would decrease remarkably if HI,were not accounted for,and scaling relation between D and c would become quite distinct.

II.MODEL AND METHOD

As shown in Fig.1(a),an NP(larger gold particle)is immersed in a polymer solution consisting of Nplinear flexible polymer chains composed of Nbbeads(green particles)of mass Mb,embedded in an explicit solvent(not shown).Exclusive volume effects among all polymer beads are taken into account by the repulsive, shifted,and truncated Lennard-Jones(LJ)potential,

where the cutoff radius is rc=21/6σ.Bond effects among adjacent beads in the polymer backbone are given by the finite extensible nonlinear elastic(FENE)potential,

where κ=30ϵ/σ2is the bond strength and R0=1.5σ is used as the maximum bond length.

The NP is modeled as an LJ-like sphere of radius Rnand mass Mn.The interaction between NP and polymer beads are described by a modified LJ potential

which is offset by the interaction range Rev=Rn−σ/2,

FIG.1(a)A typical simulation snapshot of an NP immersed in a semidilute polymer solution(with box size L=32a0). Here only nanoparticle(big yellow sphere)and polymer(red and blue beads for the both ends of a chain)is drawn,water is omitted for clarity.(b)Potential curves used in our system.Polymer beads interaction is plotted according to Eq.(1)shown with green line and nanoparticle-bead interaction is showed by orange curve referred to Eq.(3).

This potential is truncated and shifted at the separation rc=Rev+21/6σ.As an illustration,bead-bead and NP-bead interaction potentials are drawn in Fig.1(b).NP and polymer chains constitute the solutes of the system, and MD method is used to simulate their motions.

We consider that the solutes mentioned above are dissolved in the solvent,which could be water in real polymer solutions.Specific interactions among the solvent particles are not explicitly accounted for in the present work.However,we are mainly interested in the effects of HI.To this end,as already mentioned in the introduction part,we use the MPCD method which has already been widely adopted in this respect[44,47,48]. In MPCD,Nspoint solvent particles with mass ms, representing coarse-grained molecules,free stream and undergo effective collisions at discrete time intervals ΔtMPC.During streaming steps,particles move ballistically and their positions riare updated according to

where vi(t)is the velocity of particle i at time t.In collision steps,particles are sorted into cubic cells of size a0and then the velocities of solvent particles in a certain cell are updated according to a stochastic rotation

where j runs over all of the Ncparticles in the cell andis the corresponding centerof-mass velocity.is a rotation matrix which rotates velocities by an angle α around an axis generated randomly for each cell at each collision step.This simple collision rule conserves mass,momentum,and energy, which guarantees the emergence of Navier-Stokes hydrodynamics on length scales larger than the collision cell size a0and thus properly accounts for the long range HI effects.Note that a random shift of the collision lattice is necessary at every collision step to guarantee Galilean invariance[49].

In our hybrid MD-MPCD scheme,a specific coupling between the solute and solvent particles should also be taken into account.Firstly,the coupling between MPC fluid and the polymer is established in collision steps, where the polymer beads are treated similarly to solvent particles.In this case,bead size is not explicitly considered,while the momentum transfer between polymer beads and solvent particles is maintained.This kind of coupling has been widely used in polymer solutions [44,47,48].For NP,however,we may be interested in how its diffusion behavior depends on its size.To this end,the coupling between NP and solvent particles is realized by using MD simulations according to the modified LJ potential in a similar form in Eq.(3).Since the MPC solvent particle is a point particle,the offset range between solvent and NP should be Rev=Rn−σ.Note that such a similar kind of force coupling schemes has been widely used in many studies on colloidal particles immersed in a MPC solvent[50-52].

To study the effects of HI,one should investigate the same system with and without HI.A method to switch off HI in the MPC method has been proposed by Yeomans[38,39],and its basic idea is to interchange velocities of all solvent particles randomly after each collision step.Then the momentum conservation in a local field will be destroyed,and the long-ranged hydrodynamic correlations will disappear,while leaving friction coefficient and fluid self-diffusion coefficient largely unaffected.In the present work,Fisher-Yates shuffle algorithm with O(N)time complexity will be used to interchange the velocities randomly.This makes it very convenient to address the specific role of HI,by run-ning simulations with the same initial conditions and parameter values but with HI present or not.

In our simulation we set the collision cell size a0,solvent particle mass msand kBT to be the unit of length, mass,and energy,respectively.With these units,the time unit readsTo guarantee a constant local system temperature T during a simulation run,we use an efficient thermostat by rescaling the solvent velocities on the cell level as proposed in Ref.[53]. MD equations for NP and polymer beads are integrated using time-reversible velocity Verlet algorithm with a time step ΔtMD=0.002τ.The time interval between succeeding collision steps is ΔtMPC=0.1τ,the number density of solvent is ρs=10(each cell contains 10 coarsegrained particles),the mass of NP is Mn=4πR3nρsms/3, the bead size is σ=1 and the bead mass is chosen to be M=ρsms.Note that the particular settings for the bead have been shown to lead to an adequate fluidbead coupling,necessary for the appearance of hydrodynamic behavior[54].The simulation is performed in a 3-dimensional cubic box of size L=16a0.The number of monomers per chain is Nb=50 if not otherwise stated. For the accuracy of statistics,100 independent simulation runs with typically tsimu=1×108τ are performed for every system.

III.SIMULATION RESULTS AND DISCUSSION

A.Gyration radius of polymer chain

Before we study the diffusion dynamics of NP,we first investigate some relevant static properties of the polymer solution without the NP.In particular,we are interested in how the radius of gyration Rgof the polymer chain depends on the concentration c of the solution. By definition,Rgis given by

where riis the position of bead i and Rcmis the centerof-mass position of a polymer chain.We note here that there already exist well-established scaling relations regarding Rgin dilute and semi-dilute polymer solutions. For a dilute solution,there is no overlap between different polymer chains,such that the gyration radius only depends on the number of beads Nbin the chain,i.e.,

with a Flory exponent ν≈0.59 for a good solvent in theory.Note here we use a subscript‘0'to specify the value of a dilute solution which is thus not dependent on the polymer concentration.With the increase of concentration,however,polymer chains become to overlap if the monomer concentration c=NbNp/V exceeds the value given by

FIG.2(a)Log-log plot of the mean-squared radius of gyration as a function of Nb,the number of beads in a single polymer chain.The solid line represents a linear fit of the data and slope indicates ν=0.588≈0.59 according to Eq.(7). (b)Log-log plot of the relative mean-square radius of gyration as a function of the normalized concentration of the polymer solution with chain length Nb=50.The solid line represents a linear fit of the data for c＞c∗.The slope of the solid line is about-0.115 which indicates ν≈0.588 according to Eq.(9).

In the range c≫c∗,it has been shown that

for ν≈0.59.

In Fig.2(a),a log-log plot of Rg0versus Nb,obtained from our simulation for a single polymer chain in the solvent,is presented.Apparently,the scaling relation(7)is reproduced,with the Flory exponent ν≈0.58 which agrees very well with the theoretical value 0.59. In Fig.2(b),we present the relative radii of gyration Rg/Rg0as a function of the normalized concentration c/c∗.The length of the polymer chain is Nb=50 and we change Npto increase the concentration c.Clearly,for dilute solutions where c≪c∗,Rgkeeps nearly constant which is not dependent on c which is consistent with Eq.(7).For c≥c∗,a good power-law scaling relation appears between Rgand c,i.e.,Rg-c−γwith an exponent γ≈0.115 which is in excellent agreement with Eq.(9). Therefore,our simulation method can reproduce the scaling properties of the gyration radius excellently in the polymer solution,which serves as a validation of the approach in the present work.

B.NP diffusion:concentration dependence

We now turn to study the diffusion dynamics of the NP in the polymer solution.Since the number concentration of NPs in real experiments is usually very low, here we only need to consider one NP in our system and ensemble averaging is performed to get reliable data. The essential quantity is the mean square displacement (MSD)of the NP,

We can calculate the long-time diffusion coefficient D by

According to Fig.3(a),D decreases monotonically with the increment of the concentration c as expected.As proposed by Phillies[27],the dependence of D on c may follow a scaling relationship described by the stretched exponential function

where α and δ are two exponents to be determined,D0is the diffusion coefficient in a very dilute solution.He argued that this should be the case when HI dominates over topological constraints on the probe diffusion and the exponent δ range from 0.5 to 1 depending on the solvent quality.

In Fig.3(b),our simulation results of the diffusion constant D as a function of c is shown for the particle size Rn=2.The solid line is a fitting of Eq.(12)with α=4.63 and δ=0.98.Clearly,our data follows excellently the scaling relationship and the exponent δ lie within a reasonable range.In addition,we note that this result agrees very well with a recent experimental observation by Omari et al.[15]who studied the diffusion of gold NPs in semidilute solutions of polystyrene in toluene.

FIG.3(a)Mean squared distance of nanoparticles with Rn=2 for the various concentrations of polymer monomers. The short lines indicate the dependencies〈ΔR2(t)〉-t2at the beginning of time,〈ΔR2(t)〉-t for the long-time MSD, and〈ΔR2(t)〉-t0.725in the intermediate phase of high concentration case.(b)Dependence of diffusion coefficient of nanoparticle with Rn=2 on the concentration of polymer solution.The solid line is fitted via Phillies equation Eq.(12) with α=4.63,δ=0.98.

C.NP diffusion:effect of HI

As already mentioned in the Model and Method part, one of the advantages of the MPC method used here is the convenience to investigate the role of HI,by switching off HI during the simulation with all other settings the same.The basic idea is to interchange velocities of all solvent particles randomly after each collision step so that momentum(and energy)is not conserved locally. In this section,we will mainly discuss how the diffusion behavior of the NP depends on HI.

In Fig.4(a),MSD of NP as a function of time t is shown for Rn=2 and c=0.244 with HI on or off.Obviously,the long-time diffusion constant D becomes smaller(the curve is lower)if HI is switched off,indicating that HI is favorable for the NP diffusion.Interestingly,there seems to exist a sub-diffusion region in the intermediate time range when HI is not considered,and it takes longer time for the NP to enter a normal diffusion region.This indicates that HI can accelerate the relaxation of the cage formed by polymer beads surrounding NP.In Table I,the dependence of D on concentration c is listed.We note that the dif-fusion constant Don(HI on)is much larger than that Doff(HI off)in the whole concentration range considered here.In particular,the ratio Don/Doffincreases as c decreases,suggesting that the role of HI on diffusion becomes much significant in the low-concentration range.We also tried to use Eq.(12)to fit the data with HI off,which is depicted in Fig.4(b)for Rn=2.The fitting is also good,but now with exponents to be α=1.7 and δ=1.1.As discussed above,the exponent δ lies typically within the range(0.5,1)such that δ=1.1 seems not reasonable.Therefore,the effects of HI must be properly taken into account to illustrate real experimental results.

FIG.4(a)Mean squared displacement of nanoparticle with radius Rn=2 in a polymer solution with c=0.244.(b)Dependence of diffusion coefficient of nanoparticle with Rn=2 on concentration of polymer solution with HI switched off. The solid line is fitted via Phillies equation Eq.(12)with α=1.7 and δ=1.1.

TABLE I Diffusion coefficient of NP with Rn=2 as a function of polymer beads concentration,where Donand Doffrepresent the diffusion coefficients of NPs in polymer systems with HI switched on and off respectively.

FIG.5 Dependence of diffusion coefficient of NP on concentration of polymer solution with Rn=1.0,Rn=1.5,Rn=2, and Rn=2.5. Solid lines are fitted by Phillies equation Eq.(12).

TABLE II Parameters for Phillies fit in Fig.5.

D.NP diffusion:size effect

In this section,we investigate the size effects of NPs on their diffusion behaviors in polymer solutions.As already discussed before,in order to obtain reasonable mobility of NP with different size,long range HI among polymer beads has been considered in our simulation systems.Dependencies of diffusion coefficient D on polymer concentration c for a series of NPs sizes (Rn=1.0,1.5,2.0,and 2.5)are shown in Fig.5.As expected,one can see that diffusion coefficients decrease with increasing polymer concentration for all sizes of NPs considered here,but the mobility of NP slows down more remarkably for small NP than the big one.On the other hand,it can be seen that D decreases as particle size Rnincreasing for a certain polymer concentration while such a size effect gets weaker for larger polymer solution concentrations.In particular,we find that the curves can all be well-fitted by the Phillies'equation (Eq.(12)),and the fitting exponents α and δ are given in Table II.Interestingly,the exponent δ keeps nearly constant(～0.97)and does not depend on the particle size,while the exponent α increases monotonically with Rn(see Table II).We note that these simulation results are consistent with recent experimental observations in the semi-dilute concentration range[16,28].

IV.CONCLUSION

In conclusion,we used a hybrid mesoscopic MDMPCD simulation method to investigate the diffusiondynamics of NPs in semidilute polymer solutions.We demonstrated that the method can well reproduce the scaling relations between the radius of gyration and the concentration of polymer solution,which serves as a validation of the simulation approach.Extensive simulations were performed to calculate the long-time diffusion coefficient D of the NP,particularly,as functions of the solution concentration c.Interestingly,we found that the dependence of D on c can be fitt(ed ver)y well by a stretched exponential function D-exp−αcδ.The exponent α increases almost monotonically with the particle size Rn,while the exponent δ keeps a nearly constant～0.97 for the size range considered here.In addition,our method makes it convenient to investigate the very role of HI by performing simulations with HI off.It was shown that HI plays a significant role in the diffusion dynamics of the NP.On the one hand,HI can enhance the long-time diffusion constant considerably and shorten the transient time between the ballistic motion and normal diffusion.On the other hand,the scaling relation between D and c changes remarkably if HI is off and the exponent δ may become unreasonable.Our simulation results are in good agreements with recent experimental observations for NP diffusion in semidilute polymer solutions[16,28].Our results demonstrate that the hybrid MD-MPCD method is a very efficient strategy to investigate the structure and dynamics involving complex polymer solutions.

V.ACKNOWLEDGMENTS

Thiswork issupported by the NationalBasic Research Program of China(No.2013CB834606), the National Natural Science Foundation of China (No.21125313,No.21473165,No.21403204),the Foundation for Innovative Research Groups of the National Natural Science Foundation of China(No.21521001), and the Fundamental Research Funds for the Central Universities(No.WK2060030018,No.2340000034).

[1]Q.Lu and M.J.Solomon,Phys.Rev.E 66,061504 (2002).

[2]T.A.Waigh,Rep.Prog.Phys.68,685(2005).

[3]P.Cicuta and A.M.Donald,Soft Matter 3,1449 (2007).

[4]D.S.Banks and C.Fradin,Biophys.J.89,2960(2005).

[5]A.B.Goins,H.Sanabria,and M.N.Waxham,Biophys. J.95,5362(2008).

[6]T.Kalwarczyk,N.Ziebacz,A.Bielejewska,E.Zaboklicka,K.Koynov,J.Szymanski,A.Wilk,A.Patkowski, J.Gapinski,H.J.Butt,and H.Robert,Nano Lett.11, 2157(2011).

[7]M.Tabaka,T.Kalwarczyk,J.Szymanski,S.Hou,and R.Holyst,Biophysics 2,54(2014).

[8]J.Liu,L.Zhang,D.Cao,and W.Wang,Phys.Chem. Chem.Phys.11,11365(2009).

[9]S.K.Kumar and R.Krishnamoorti,Annual Rev. Chem.Biomol.Eng.1,37(2010).

[10]B.Hong and A.Z.Panagiotopoulos,Soft Matter 9, 6091(2013).

[11]Y.Li,M.Kr¨oger,and W.K.Liu,Soft Matter 10,1723 (2014).

[12]A.Karatrantos,N.Clarke,R.J.Composto,and K.I. Winey,Soft Matter 11,382(2015).

[13]X.Ye,P.Tong,and L.Fetters,Macromolecules 31, 5785(1998).

[14]R.Holyst,A.Bielejewska,J.Szyma´nski,A.Wilk,A. Patkowski,J.Gapi’nski,A.BZywoci´nski,T.Kalwarczyk,E.Kalwarczyk,M.Tabaka,Z.Natalia,and W. A.Stefan,Phys.Chem.Chem.Phys.11,9025(2009).

[15]R.A.Omari,A.M.Aneese,C.A.Grabowski,and A. Mukhopadhyay,J.Phys.Chem.B 113,8449(2009).

[16]I.Kohli and A.Mukhopadhyay,Macromolecules 45, 6143(2012).

[17]H.Guo,G.Bourret,R.B.Lennox,M.Sutton,J.L. Harden,and R.L.Leheny,Phys.Rev.Lett.109,055901 (2012).

[18]E.A.Mun,C.Hannell,S.E.Rogers,P.Hole,A.C. Williams,and V.V.Khutoryanskiy,Langmuir 30,308 (2013).

[19]A.Vagias,R.Raccis,K.Koynov,U.Jonas,H.J.Butt, G.Fytas,P.Ko˘sovan,O.Lenz,and C.Holm,Phys. Rev.Lett.111,088301(2013).

[20]F.Babaye Khorasani,R.Poling-Skutvik,R.Krishnamoorti,and J.C.Conrad,Macromolecules 47,5328 (2014).

[21]R.Poling-Skutvik,R.Krishnamoorti,and J.C.Conrad, ACS Macro Lett.4,1169(2015).

[22]L.H.Cai,S.Panyukov,and M.Rubinstein,Macromolecules 44,7853(2011).

[23]A.Ochab-Marcinek and R.Holyst,Soft Matter 7,7366 (2011).

[24]A.Ochab-Marcinek,S.A.Wieczorek,N.Ziebacz,and R.Holyst,Soft Matter 8,11173(2012).

[25]U.Yamamoto and K.S.Schweizer,Macromolecules (2014).

[26]L.H.Cai,S.Panyukov,and M.Rubinstein,Macromolecules 48,847(2015).

[27]G.D.Phillies,Macromolecules 20,558(1987).

[28]A.Michelman-Ribeiro,F.Horkay,R.Nossal,and H. Boukari,Biomacromolecules 8,1595(2007).

[29]R.Tuinier and T.H.Fan,Soft Matter 4,254(2008).

[30]S.Egorov,J.Chem.Phys.134,084903(2011).

[31]U.Yamamoto and K.S.Schweizer,J.Chem.Phys.135, 224902(2011).

[32]X.Y.Lai and N.R.Zhao,Chin.J.Chem.Phys.26, 163(2013).

[33]Y.Dong,X.Feng,N.Zhao,and Z.Hou,J.Chem.Phys. 143,024903(2015).

[34]J.Liu,Y.Gao,D.Cao,L.Zhang,and Z.Guo,Langmuir 27,7926(2011).

[35]J.Liu,D.Cao,and L.Zhang,J.Phys.Chem.C 112, 6653(2008).

[36]T.K.Patra and J.K.Singh,J.Chem.Phys.138, 144901(2013).

[37]J.T.Kalathi,U.Yamamoto,K.S.Schweizer,G.S. Grest,and S.K.Kumar,Phys.Rev.Lett.112,108301 (2014).

[38]N.Kikuchi,A.Gent,and J.Yeomans,European Phys.J.E 9,63(2002).

[39]N.Kikuchi,J.Ryder,C.Pooley,and J.Yeomans,Phys. Rev.E 71,061804(2005).

[40]R.Chang and A.Yethiraj,J.Chem.Phys.114,7688 (2001).

[41]H.Tanaka,J.Phys.:Condens.Matter 17,S2795 (2005).

[42]K.Kamata,T.Araki,and H.Tanaka,Phys.Rev.Lett. 102,108303(2009).

[43]A.Malevanets and R.Kapral,J.Chem.Phys.110,8605 (1999).

[44]A.Malevanets and J.Yeomans,EPL 52,231(2000).

[45]R.Kapral,Adv.Chem.Phys.140,89(2008).

[46]G.Gompper,T.Ihle,D.Kroll,and R.Winkler,in Advanced Computer Simulation Approaches for Soft Matter Sciences III,Berlin,Heidelberg:Springer,1(2009).

[47]C.C.Huang,R.G.Winkler,G.Sutmann,and G. Gompper,Macromolecules 43,10107(2010).

[48]C.Echeverria and R.Kapral,Phys.Chem.Chem.Phys. 14,6755(2012).

[49]T.Ihle and D.M.Kroll,Phys.Rev.E 63,020201 (2001).

[50]A.Malevanets and R.Kapral,J.Chem.Phys.112,7260 (2000).

[51]J.Padding and A.Louis,Phys.Rev.Lett.93,220601 (2004).

[52]A.Tomilov,A.Videcoq,T.Chartier,T.Ala-Nissil¨a, and I.Vattulainen,J.Chem.Phys.137,014503(2012).

[53]C.Huang,A.Chatterji,G.Sutmann,G.Gompper,and R.G.Winkler,J.Comput.Phys.229,168(2010).

[54]M.Ripoll,K.Mussawisade,R.G.Winkler,and G. Gompper,Phys.Rev.E 72,016701(2005).

(Dated:Received on March 27,2016;Accepted on May 4,2016)