The article information
 Li Shuxian, Jiang Huijun, Hou Zhonghuai
 李树贤, 江慧军, 侯中怀
 Diffusion of Nanoparticles in Semidilute Polymer Solutions: A Multiparticle Collision Dynamics Study
 纳米粒子在亚浓高分子溶液中的扩散: 多粒子碰撞动力学
 Chinese Journal of Chemical Physics, 2016, 29(5): 549556
 化学物理学报, 2016, 29(5): 549556
 http://dx.doi.org/10.1063/16740068/29/cjcp1603058

Article history
 Received on: March 27, 2016
 Accepted on: May 4, 2016
b. Department of Chemical Physics and Hefei National Laboratory for Physical Sciences at the Microscales,
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 ^{[13]}. 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 ^{[47]}. 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 ^{[812]}, to list just a few.
For many of these reasons, diffusion of NPs in polymer solutions has received a lot of attention both experimentally ^{[1321]} and theoretically ^{[2226]}. It was found that the wellknown StokesEinstein (SE) relationship between the longtime 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 semidilute transition regime, the dependence of D on c was shown to follow a stretched exponential function given by the Phillies equation ^{[27]}, Dexp (−α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, 3033]}, 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 longtime diffusion coefficient in complex polymer liquids, from a microscopic level explicitly accounting for the solutesolvent 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 longtime transverse hydrodynamic currents. With the calculated ζ(t), the diffusion coefficient D can be obtained directly via fluctuationdissipation 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 coarsegrained molecular dynamics simulations (CGMD) as well ^{[11, 3437]}. 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]}.
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 longrange hydrodynamic interaction (HI), might drastically influence dynamics of polymer chain in a solvent ^{[3842]}. 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 coarsegrained 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 semidilute 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 semidilute 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 METHODAs shown in Fig. 1(a), an NP (larger gold particle) is immersed in a polymer solution consisting of N_{p} linear flexible polymer chains composed of N_{b} beads (green particles) of mass M_{b} , embedded in an explicit solvent (not shown). Exclusive volume effects among all polymer beads are taken into account by the repulsive, shifted, and truncated LennardJones (LJ) potential,
(1) 
where the cutoff radius is r_{c}=2^{1/6}σ. Bond effects among adjacent beads in the polymer backbone are given by the finite extensible nonlinear elastic (FENE) potential,
(2) 
where κ=30ϵ/σ^{2} is the bond strength and R_{0}=1. 5σ is used as the maximum bond length.
The NP is modeled as an LJlike sphere of radius R_{n} and mass Mn. The interaction between NP and polymer beads are described by a modified LJ potential which is offset by the interaction range R_{ev}=R_{n}−σ/2,
(3) 
This potential is truncated and shifted at the separation r_{c}=R_{ev}+2^{1/6}σ. As an illustration, beadbead and NPbead 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, N_{s} point solvent particles with mass m_{s}, representing coarsegrained molecules, free stream and undergo effective collisions at discrete time intervals Δt_{MPC}. During streaming steps, particles move ballistically and their positions r_{i} are updated according to
(4) 
where v_{i}(t) is the velocity of particle i at time t. In collision steps, particles are sorted into cubic cells of size a_{0} and then the velocities of solvent particles in a certain cell are updated according to a stochastic rotation
(5) 
where j runs over all of the Nc particles in the cell and
In our hybrid MDMPCD 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 R_{ev}=R_{n}−σ. 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 ^{[5052]}.
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 longranged hydrodynamic correlations will disappear, while leaving friction coefficient and fluid selfdiffusion coefficient largely unaffected. In the present work, FisherYates 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 running simulations with the same initial conditions and parameter values but with HI present or not.
In our simulation we set the collision cell size a_{0}, solvent
particle mass m_{s} and k_{B}T to be the unit of length,
mass, and energy, respectively. With these units, the
time unit reads
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 R_{g} of the polymer chain depends on the concentration c of the solution. By definition, R_{g} is given by
(6) 
where r_{i} is the position of bead i and Rcm is the centerofmass position of a polymer chain. We note here that there already exist wellestablished scaling relations regarding R_{g} in dilute and semidilute 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 N_{b} in the chain, i. e. ,
(7) 
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=N_{b}N_{p}/V exceeds the value given by
(8) 
In the range c≫c^{*}, it has been shown that
(9) 
for ν≈0. 59. In Fig. 2(a), a loglog plot of R_{g}0 versus N_{b}, 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 R_{g}/R_{g}0 as a function of the normalized concentration c/c^{*}. The length of the polymer chain is N_{b}=50 and we change N_{p} to increase the concentration c. Clearly, for dilute solutions where c≪c^{*}, R_{g} keeps nearly constant which is not dependent on c which is consistent with Eq. (7). For c≥c^{*}, a good powerlaw scaling relation appears between R_{g} and c, i. e. , R_{g}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 dependenceWe 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,
(10) 
where r_{n} (t) denotes the position vector of NP at time t and
We can calculate the longtime diffusion coefficient D by
(11) 
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
(12) 
where α and δ are two exponents to be determined, D0 is 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 R_{n}=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.
C. NP diffusion: effect of HIAs 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 R_{n}=2 and c=0. 244 with HI on or off. Obviously, the longtime 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 subdiffusion 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 diffusion constant Don (HI on) is much larger than that Doff (HI off) in the whole concentration range considered here. In particular, the ratio Don/Doff increases as c decreases, suggesting that the role of HI on diffusion becomes much significant in the lowconcentration range. We also tried to use Eq. (12) to fit the data with HI off, which is depicted in Fig. 4(b) for R_{n}=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.
D. NP diffusion: size effectIn 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 (R_{n}=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 R_{n} increasing 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 wellfitted 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 R_{n} (see Table II). We note that these simulation results are consistent with recent experimental observations in the semidilute concentration range ^{[16, 28]}.
IV. CONCLUSIONIn conclusion, we used a hybrid mesoscopic MDMPCD simulation method to investigate the diffusion dynamics 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 longtime 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 fitted very well by a stretched exponential function Dexp (−αc^{δ} ) . The exponent α increases almost monotonically with the particle size R_{n}, 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 longtime 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 MDMPCD method is a very efficient strategy to investigate the structure and dynamics involving complex polymer solutions.
V. ACKNOWLEDGMENTSThis work is supported by the National Basic 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]  Lu M. Q.,and Solomon J., Phys. Rev E66 , 061504 (2002). 
[2]  A. Waigh T., Rep. Prog. Phys. 68 , 685 (2005). DOI:10.1088/00344885/68/3/R04 
[3]  Cicuta A. P.,and Donald M., Soft Matter 3 , 1449 (2007). DOI:10.1039/b706004c 
[4]  S. Banks D.,and Fradin C., Biophys. J. 89 , 2960 (2005). DOI:10.1529/biophysj.104.051078 
[5]  B. Goins A., Sanabria H.,and N.Waxham M., Biophys. J. 95 , 5362 (2008). DOI:10.1529/biophysj.108.131250 
[6]  Kalwarczyk T., Ziebacz N., Bielejewska A., Zaboklicka E., Koynov K., Szymanski J., Wilk A., Patkowski A., Gapinski J., J. Butt H.,and Robert H., Nano Lett. 11 , 2157 (2011). DOI:10.1021/nl2008218 
[7]  Tabaka M., Kalwarczyk T., Szymanski J., Hou S.,and Ho lyst R., Biophysics 2 , 54 (2014). 
[8]  Liu J., Zhang L., Cao D.,and Wang W., Phys. Chem. Chem. Phys. 11 , 11365 (2009). DOI:10.1039/b913511a 
[9]  K. Kumar R. Krishnamoorti S., Annual Rev. Chem. Biomol. Eng. 1 , 37 (2010). DOI:10.1146/annurevchembioeng073009100856 
[10]  Hong A. Z. Panagiotopoulos B., Soft Matter 9 , 6091 (2013). DOI:10.1039/c3sm50832c 
[11]  Li Y., Kröger M.,and K. Liu W., Soft Matter 10 , 1723 (2014). DOI:10.1039/c3sm51564h 
[12]  Karatrantos A., Clarke N., J. Composto R.,and I. Winey K., Soft Matter 11 , 382 (2015). DOI:10.1039/C4SM01980F 
[13]  Ye X., Tong P.,and Fetters L., Macromolecules 31 , 5785 (1998). DOI:10.1021/ma9801725 
[14]  Ho lyst R., Bielejewska A., Szymański J., Wilk A., Patkowski A., Gapi'nski J., BZywociński A., Kalwarczyk T., Kalwarczyk E., Tabaka M., Natalia Z.,and A. Stefan W., Phys. Chem. Chem. Phys. 11 , 9025 (2009). DOI:10.1039/b908386c 
[15]  A. Omari R., M. Aneese A., A. Grabowski C.,and Mukhopadhyay A., J. Phys. Chem B113 , 8449 (2009). 
[16]  Kohli A. Mukhopadhyay I., Macromolecules 45 , 6143 (2012). DOI:10.1021/ma301237r 
[17]  Guo H., Bourret G., B. Lennox R., Sutton M., L. Harden J.,and L. Leheny R., Phys. Rev. Lett. 109 , 055901 (2012). DOI:10.1103/PhysRevLett.109.055901 
[18]  A. Mun E., Hannell C., E. Rogers S., Hole P., C. Williams A.,and V. Khutoryanskiy V., Langmuir 30 , 308 (2013). 
[19]  Vagias A., Raccis R., Koynov K., Jonas U., J. Butt H., Fytas G., Košovan P., Lenz O.,and Holm C., Phys. Rev. Lett. 111 , 088301 (2013). DOI:10.1103/PhysRevLett.111.088301 
[20]  Babaye Khorasani F., PolingSkutvik R., Krishnamoorti R.,and C. Conrad J., Macromolecules 47 , 5328 (2014). DOI:10.1021/ma501248u 
[21]  PolingSkutvik R., Krishnamoorti R.,and C. Conrad J., ACS Macro Lett. 4 , 1169 (2015). DOI:10.1021/acsmacrolett.5b00616 
[22]  H. Cai L., Panyukov S.,and Rubinstein M., Macromolecules 44 , 7853 (2011). DOI:10.1021/ma201583q 
[23]  OchabMarcinek R. Ho lyst A., Soft Matter 7 , 7366 (2011). DOI:10.1039/c1sm05217a 
[24]  OchabMarcinek A., A. Wieczorek S., Ziebacz N.,and Ho lyst R., Soft Matter 8 , 11173 (2012). DOI:10.1039/c2sm25925g 
[25]  U. Yamamoto and K. S. Schweizer, Macromolecules (2014). 
[26]  H. Cai L., Panyukov S.,and Rubinstein M., Macromolecules 48 , 847 (2015). DOI:10.1021/ma501608x 
[27]  D. Phillies G., Macromolecules 20 , 558 (1987). DOI:10.1021/ma00169a015 
[28]  MichelmanRibeiro A., Horkay F., Nossal R.,and Boukari H., Biomacromolecules 8 , 1595 (2007). DOI:10.1021/bm061195r 
[29]  Tuinier T. H. Fan R., Soft Matter 4 , 254 (2008). DOI:10.1039/B711902J 
[30]  Egorov S., J. Chem. Phys. 134 , 084903 (2011). DOI:10.1063/1.3556749 
[31]  Yamamoto K. S. Schweizer U., J. Chem. Phys. 135 , 224902 (2011). DOI:10.1063/1.3664863 
[32]  Y. Lai N. R. Zhao X., Chin. J. Chem. Phys. 26 , 163 (2013). DOI:10.1063/16740068/26/02/163171 
[33]  Dong Y., Feng X., Zhao N.,and Hou Z., J. Chem. Phys. 143 , 024903 (2015). DOI:10.1063/1.4926412 
[34]  Liu J., Gao Y., Cao D., Zhang L.,and Guo Z., Langmuir 27 , 7926 (2011). DOI:10.1021/la201073m 
[35]  Liu J., Cao D.,and Zhang L., J. Phys. Chem C112 , 6653 (2008). 
[36]  K. Patra J. K. Singh T., J. Chem. Phys. 138 , 144901 (2013). DOI:10.1063/1.4799265 
[37]  T. Kalathi J., Yamamoto U., S. Schweizer K., S. Grest G.,and K. Kumar S., Phys. Rev. Lett. 112 , 108301 (2014). DOI:10.1103/PhysRevLett.112.108301 
[38]  Kikuchi N., Gent A.,and Yeomans J., European Phys. J E9 , 63 (2002). 
[39]  Kikuchi N., Ryder J., Pooley C.,and Yeomans J., Phys. Rev E71 , 061804 (2005). 
[40]  Chang A. Yethiraj R., J. Chem. Phys. 114 , 7688 (2001). DOI:10.1063/1.1361071 
[41]  Tanaka H., J. Phys.: Condens. Matter 17 , S2795 (2005). DOI:10.1088/09538984/17/31/004 
[42]  Kamata K., Araki T.,and Tanaka H., Phys. Rev. Lett. 102 , 108303 (2009). DOI:10.1103/PhysRevLett.102.108303 
[43]  Malevanets R. Kapral A., J. Chem. Phys. 110 , 8605 (1999). DOI:10.1063/1.478857 
[44]  Malevanets J. Yeomans A., EPL 52 , 231 (2000). DOI:10.1209/epl/i2000004280 
[45]  Kapral R., Adv. Chem. Phys. 140 , 89 (2008). DOI:10.1002/SERIES2007 
[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. Huang C., G. Winkler R., Sutmann G.,and Gompper G., Macromolecules 43 , 10107 (2010). DOI:10.1021/ma101836x 
[48]  Echeverria R. Kapral C., Phys. Chem. Chem. Phys. 14 , 6755 (2012). DOI:10.1039/c2cp40200a 
[49]  Ihle D. M. Kroll T., Phys. Rev E63 , 020201 (2001). 
[50]  Malevanets R. Kapral A., J. Chem. Phys. 112 , 7260 (2000). DOI:10.1063/1.481289 
[51]  Padding A. Louis J., Phys. Rev. Lett. 93 , 220601 (2004). DOI:10.1103/PhysRevLett.93.220601 
[52]  Tomilov A., Videcoq A., Chartier T., AlaNissilä T.,and Vattulainen I., J. Chem. Phys. 137 , 014503 (2012). DOI:10.1063/1.4731661 
[53]  Huang C., Chatterji A., Sutmann G., Gompper G.,and G. Winkler R., J. Comput. Phys. 229 , 168 (2010). DOI:10.1016/j.jcp.2009.09.024 
[54]  Ripoll M., Mussawisade K., G. Winkler R.,and Gompper G., Phys. Rev E72 , 016701 (2005). 
b. 中国科学技术大学合肥微尺度物质科学国家实验室(筹), 合肥 230026