Dissipative Particle Dynamics Simulations of Grafted Polymer Chains

4732

Langmuir 2000, 16, 4732-4740

Dissipative Particle Dynamics Simulations of Grafted Polymer Chains between Two Walls P. Malfreyt Laboratoire de thermodynamique des Solutions et des Polyme` res UPRES A 6003, 24 avenue des Landais, Universite´ Blaise-Pascal, 63177 Aubie` re Cedex, France

D. J. Tildesley* Unilever Research Port Sunlight, Bebington, Wirral L63 3JW, U.K. Received October 22, 1999. In Final Form: January 10, 2000 We have used the dissipative particle dynamics method (DPD) to simulate grafted polymer brushes in a pore geometry using two walls of tethered DPD particles to form the top and bottom boundaries of the simulation cell. The polymer density profiles show a parabolic shape in agreement with the SCF theory for brushes in a good solvent. At fixed wall separation and chain length, the overlap of the brushes increases with increasing surface coverage. The widths of the polymer brushes decrease significantly with decreasing surface coverage, and the effects are more dramatic for the longer chains. The diffusion along the pore axis is significantly greater for the solvent particles in the middle of the slab than for those in the polymer/wall region since the solvent molecules are trapped in the entangled polymer. Calculation of the average asphericity parameters from the inertia tensor of the polymers demonstrates that the chains stretch with increasing coverage. As the chains become longer, the degree of ordering with respect to the layer normal increases at a given surface density. The ordering with respect to the surface normal decreases with decreasing coverage at fixed chain length. The director of the layer is close to the layer normal.

I. Introduction Polymers grafted or adsorbed on a surface give rise to a wide range of important industrial applications in automotive lubrication, oil recovery, modification of surfaces, friction, and contact formation. In particular understanding the lubrication process requires a detailed study of the effect of these grafted polymer brushes on the rheology of a liquid in a confined geometry. In this paper, we present a simulation study of grafted brushes in a good solvent confined between fixed walls. The relaxation of these brushes occurs on the time scales longer than 10-9 s, which makes the use of the standard molecular dynamics method with an atomistic model inappropriate. The principal purpose of this paper is to develop a coarsegrained model of the system and to apply a new mesoscopic simulation technique, dissipative particle dynamics (DPD), to investigate the equilibrium behavior of the brushes. This is a necessary prerequisite to a more detailed modeling study of the shearing of polymer brushes and a direct calculation of the friction coefficient and viscosity. There have been a significant number of experimental and theoretical studies on the properties of polymer brushes. Early theories by Alexander1 and de Gennes2 suggested that at high grafting density the polymer chains are strongly stretched in the direction normal to the grafting surface. In an extension of the work of Alexander, Milner et al.3 demonstrated that a balance between the chain extension and the bead interaction resulted in a parabolic density profile for a single chain in the mean field approximation. The parabolic profile is consistent with experimental results4 and the results of simula(1) Alexander, S. J. Phys. (Paris) 1977, 38, 983. (2) de Gennes, P. G. Scaling concepts in polymer physics; Cornell University Press: Ithaca, NY, 1979; Macromolecules 1980, 13, 1069. (3) Milner, S. T.; Witten, T. A.; Cates, M. E. Macromolecules 1988, 21, 2610. (4) Patel, S. S.; Tirrell, M. Annu. Rev. Phys. Chem. 1989, 40, 597.

tion.5-10 Grest has shown that at high grafting density, the brush height scales as FA1/2 rather than FA1/3 where FA is the grafting density of polymers on the surface. Much of the reported simulation work in this area has been carried using molecular dynamics.11,12 The method normally uses pair additive forces to describe the atomic interactions between and within molecules. The coupled equations of motion are solved numerically using a finite difference method with a time step ranging from 10-16 to 10-14 s. This time step is chosen to track the fastest motion in the system and needs to be in this range to ensure conservation energy in the simulation. This places an effective upper limit of approximately 20 ns on the real time that can be routinely simulated in this way. In many cases, the inclusion of full atomistic detail in simulations of a polymer system is unnecessary. Many physical properties of polymers can be successfully modeled using a more coarse-grained description. This description has already been successfully applied in Monte Carlo13-18 and (5) Grest, G. S.; Murat, M. In Monte Carlo and molecular dynamics simulation in polymer science; Binder, K., Ed.; Oxford University Press: London, 1995; Chapter 9. (6) Lai, P. Y.; Zhulina, E. B. Macromolecules 1992, 25, 5201. (7) Dickman, R.; Anderson, P. E. J. Chem. Phys. 1993, 99, 3112. (8) Lai, P. Y.; Binder, K. J. Chem. Phys. 1993, 98, 2366. (9) Murat, M.; Grest, G. S. Phys. Rev. Lett. 1989, 63, 1074. (10) Murat, M.; Grest, G. S. In Computer Simulations of Polymers; Roe, R. J., Ed.; Prentice-Hall: Englewood Cliffs, NJ, 1991. (11) Allen, M. P.; Tildesley, D. J. Computer simulation of Liquids; Clarendon: Oxford, U.K., 1987. (12) Allen, M. P.; Tildesley, D. J. Computer simulations in Chemical Physics; Kluwer Academic: Dordrecht, The Netherlands, 1993. (13) Baumga¨rtner, A. Application of the Monte Carlo Method in Statistical Physics; Binder, K., Ed.; Springer: Berlin, 1984. (14) Binder, K. Computational Modeling of Polymers; J. Biorano: New York, 1992. (15) Binder, K. Monte Carlo Methods in Statistical Physics; Springer: Berlin, 1979. (16) Diestler, D. J.; Schoen, M.; Cushman, J. H. Science 1993, 262, 545.

10.1021/la991396z CCC: $19.00 © 2000 American Chemical Society Published on Web 04/21/2000

Dissipative Particle Dynamics Simulations

Langmuir, Vol. 16, No. 10, 2000 4733

molecular dynamics simulations19-25 of polymers by using bead that represent the persistence or Kuhn length of the polymer and by attaching the segment using simple springs. This model can be simulated using the molecular dynamics method with a larger time step. The reduced time step in the simulation is scaled by (mL2/)1/2 where L and are the length and energy parameters associated with a segment in the bead-spring model. In this paper, we wish to explore the application of a novel mesoscale method to study lubrication. The DPD method was introduced by Hoogerbrugge and Koelman26,27 and extended to polymers by Schlijper.28,29 The method uses a soft conservative force to represent the volume element of the polymer or solvent. The repulsion is supplemented by pair-additive, Brownian, and hydrodynamic interactions that act as a thermostat to keep the system at constant temperature. Espano˜l and Warren30 have shown that the method satisfies detailed balance and at equilibrium simulates a system in the canonical ensemble. The method also conserves Newton’s third law and the dynamics exhibits Navier Stokes behavior. Typically the distance scale in the DPD method is of the order of 100 nm and the time step of the order of 25 ns. The length and time step in the method make it an attractive method to study the equilibrium structure and rheology of grafted polymer chains. However, the method relies on the softness of the conservative force to achieve this advantage. (Previous molecular dynamics simulations of the bead-spring model have used a shifted and truncated Lennard-Jones potential to represent the bead with an inverse 12th power repulsion). It is important to establish that the DPD method can reproduce the important features of the equilibrium polymer brush behavior before we can use it with confidence to study the rheology. In section 2, we outline the model and simulation. In section 3 work, we calculate the temperature and pressure profiles across the pore and report on the translational, orientational and conformational structure of the brushes as a function of surface coverage and polymer chain length. Section 4 contains our conclusions.

the chain that are equivalent to the persistence length of the polymer; the solvent particle represents an element of the fluid containing hundreds of solvent molecules. These particles are best though of as momentum carriers and the relationship between the DPD repulsion parameters, aij, and the Flory-Huggins χ parameters has been discussed by Groot and Warren31 The dissipative particle dynamics method solves the Newtonian equations of motions for particles subject to conservative, dissipative and random force. Thus

II. The Dissipative Particle Dynamics Model and Simulation

dri dvi ) vi, mi ) fi dt dt

The model used in the simulations described in this paper consists of solvent particles, polymer beads and two boundary walls. The initial geometry is sketched in Figure 1a. The DPD “particle” is a coarse-grained representation of the underlying atomistic model. The beads of the polymer represent the number of atomic repeat units of (17) Atkinson, J.; Goth, C. J.; Phan-Thien, N. J. Chem. Phys. 1984, 80, 6305. (18) Lai, P. Y.; Binder, K. J. Chem. Phys. 1993, 98, 2366. (19) Grest, G. S.; Kremer, K. Phys. Rev. 1986, A33, 3628. (20) Kremer. K.; Grest, G. S. J. Chem. Phys. 1990, 92, 5057. (21) Kremer, K.; Grest G. S In Computer Simulations of Polymers; Roe, R. J., Ed.; Prentice-Hall: Englewood, Cliffs, NJ, 1991. (22) Kremer, K. In Computer Simulation in Chemical Physics, Vol. 397 of NATO Advanced Study Institute, Series C: Mathematical And Physical Science; Allen M. P., Tildesley, D. J., Eds; Kluwer Academic: London, 1993. (23) Murat, M.; Grest, G. S. Phys. Rev. Lett. 1989, 63, 1074. (24) Peters, G. H.; Tildesley, D. J. Phys. Rev. E 1995, 52, 1882. (25) Peters, G. H.; Tildesley, D. J. Phys. Rev. E 1996, 54, 5493. (26) Hoogerbrugge, P. J.; Koelman, J. M. V. A Europhys. Lett. 1992, 19, 155. (27) Koleman, J. M. V. A.; Hoogerbrugge, P. J. Europhys. Lett. 1993, 21, 363. (28) Schlijper, A. G.; Hoggerbrugge, P. J.; Manke, C. W. J. Rheol. 1995, 39, 567. (29) Madden, W. G.; Kong, Y.; Manke; C. M.; Schlijper, A. G. Int. J. Thermophys. 1994, 15, 1093. (30) Espano˜l, P.; Warren, P. B. Europhys. Lett. 1995, 30, 191.

Figure 1. (a) Schematic representation of the simulation box. Open circles represent wall particles, circles marked with a tilde represent solvent particles, and black circles represent polymer beads. (b) Part of the first layer of wall particles viewed from above. The open circles represent wall particles that are uncovered. The black circles represent the heads of polymer chains, and the arrows represent the direction of the first polymer bond in the initial configuration.

(1)

where fi is a pairwise-additive force

fi )

R (FCij + FD ∑ ij + Fij ) j*i

(2)

FCij is a slowly varying conservative force acting between particle i and j representing the particle-particle interaction and the spring interaction sustaining the polymer R topology. FD ij and Fij are the dissipative force and random force, respectively. All interactions are calculated with a cutoff radius rc ) 1. All the distances are reduced with respect to rc and all particles have equal mass, m ) 1. The interparticle, conservative force is linear

FCij )

{

aij (1 - rij) rˆ ij 0

(rij < 1) (rij g 1)

(3)

where aij is the, reduced repulsion parameter and rij ) ri - rj, rˆ ij ) rij/|rij|. The repulsive parameter aij is set to 60 throughout the paper. The solvent particles see one another, the wall particles and the polymer particles with the same interaction. Thus, the polymer chains are (31) Groot, R. D.; Warren, P. B. J. Chem. Phys 1997, 107, 4423.

4734

Langmuir, Vol. 16, No. 10, 2000

Malfreyt and Tildesley

unlikely to adsorb onto the wall because of entropic forces and the fact that the polymer is in a good solvent. The dissipative force depends on the relative particle velocities vij ) vi - vj and is D FD ˆ ij‚vij)rˆ ij ij ) -γω (rij)(r

(4)

The random force is

FRij ) σωR(rij)θij

1 rˆ ij x∆t

(5)

In these equations γ can be interpreted as the friction coefficient and σ the amplitude of the noise. θij is a random number sampled from a Gaussian distribution with zero mean and unit variance and is chosen independently for each pair of particles at each time step. Espan˜ol and Warren30 have shown that the system will sample the canonical ensemble and obey the fluctuation-dissipation theorem if

γ)

σ2 and ωD(rij) ) [ωR(rij)]2 2kBT

(6)

Where appropriate the DPD particles are connected to form a polymer using the conservative spring force

FSij ) kS (|ri - rj| - req)rˆ ij

(7)

where kS is the spring constant and req the equilibrium spring length. The spring constant is adjusted so that the average spring length is equal to the mean nearestneighbor distance between solvent particles. The values kS ) -225 and req ) 0.85 produced an average intrapolymer distance that coincides approximately with the first maximum of the two-dimensional radial distribution function for solvent particles. The wall particles are initially arranged on a face-centered cubic lattice with the 001 surface exposed. Each wall consists of 486 atoms arranged in three layers. The nearest-neighbor separation in the lattice is 0.583 and surface consists of 9 × 9 units cells each containing two particles. The wall particles are attached to their initial lattice sites, req,i, by a force

ˆ ieq,i Fw i ) kw(|ri - req,i|)r

(8)

where kw ) -150. This value of the spring constant is sufficient to keep the wall particles in a solid structure. This solid persists throughout the simulation and prevents the polymer chains from penetrating the wall significantly. The simulations are performed in a parallelepiped box with periodic boundary conditions applied in the x and y directions; in the z direction the fluid is confined by two walls. The reduced cell dimensions are Lz ) 18.78, Lx ) 7.42, and Ly ) 7.42. Simulations are performed at coverages of 1/3 or 1/6, with 27 or 54 polymer chains respectively attached to each wall. The polymers are attached by a spring to a particular wall particle. The chains are grafted uniformly to the surface as shown in Figure 1b. The first particle of the polymer chain is part of the wall. The polymer chain length is varied between 5 and 20 beads. The total number of particles in the box is fixed at 5006 which gives an overall reduced density, Frc3 ) 4.83. Particles that do not belong to the polymer chains or the wall are considered to be solvent particles. Thus, the total density is unaffected by the chain length. These solvent particles are initially placed on an fcc lattice constructed in the middle of the box, avoiding direct

Table 1. Details of the DPD Simulations Presented in the Paper, Where the Averages Presented in This Table Are Taken over the Final 20 000 Steps of the Run length of layer chain, width, Nb coverage w ϑdirector 〈P1,n〉 〈P2,d〉 5 10 15 20 5 10 15 20

1/3 1/3 1/3 1/3 1/6 1/6 1/6 1/6

3.08 5.59 7.91 9.72 2.95 5.13 6.85 8.74

0.0949 0.0733 0.0518 0.0485 0.171 0.1192 0.1166 0.0908

0.736 0.790 0.827 0.852 0.691 0.735 0.759 0.803

0.389 0.500 0.580 0.632 0.310 0.391 0.445 0.531

〈I(LZ)〉

〈A〉

0.0000 0.0000 0.00026 0.013 0.0000 0.0000 0.000030 0.0019

0.520 0.547 0.588 0.627 0.501 0.484 0.514 0.532

overlap with the polymer chain. One advantage of the DPD method is that the soft conservative force means that particles cannot overlap strongly in any configuration and that any starting configurations can be used without prejudice. All traces of this initial configuration of solvent particles are lost in the early stages of the run. The method used to integrate the equations of motion is a modified version of the velocity-Verlet algorithm.11

1 ri(t + ∆t) ) ri(t) + ∆t vi(t) + (∆t)2fi(t) 2 1 vi(t + 1/2∆t) ) vi(t) + (∆t) fi(t) 2

(9)

1

fi(t + ∆t) ) fi(ri(t + ∆t), vi(t + /2∆t)) 1 vi(t + ∆t) ) vi(t + 1/2∆t) + (∆t) fi(t + ∆t). 2 A time step of ∆t ) 0.02 is used throughout. A Verlet neighbor list is implemented in the code using a list radius r1 ) 2.6; the neighbor list is updated automatically when the sum of the magnitude of the two largest displacements since the previous update exceeds r1 - rc. In the simulations the amplitude of the noise, σ, was chosen so that σ/x∆t ) 15 with a time step equal to 0.02. Hence, the value of the friction coefficient γ ) 1.12 corresponds to a reduced temperature T ) 2.0. The average temperature of each wall calculated from the MaxwellBoltzmann distribution of velocity components in the same as the temperature fixed from the ratio of σ and γ in eq 6. A typical simulation run consisted of 100 000 steps of equilibration followed by a production phase of additional 100 000 steps. The thermodynamic properties (such as pRβ(z), E, T (z)) and the structural quantities (such as F(z), g(s) and the radii of gyration) were calculated directly during the course of the simulation. III. Results and Discussions Table 1 contains a summary of the simulations discussed in this paper. The average temperature profiles, Tx(z), Ty(z), Tz(z), were calculated to ensure that the systems were in thermal equilibrium

〈 〉 N

TR(z) )

Hn(zi) vRi‚vRi ∑ i)1 N

(10)

Hn(zi) ∑ i)1

where 〈 〉 denotes a time average. vRi is the R-component


Langmuir, Vol. 16, No. 10, 2000 4735

Figure 2. Temperature profiles of the three components Tx(z), Ty(z), Tz(z) at a coverage of 1/3 for a polymer chain length of 15 beads. The individual temperature components are offset by 0.5 temperature units for clarity.

of the velocity of particle i, and Hn(zi) is a top-hat function

{

∆z ∆z < zi < zn + 2 2 Hn(zi) ) 0 otherwise Hn(zi) ) 1 for zn -

(11)

The function was evaluated at intervals of 10 time steps throughout the production phase of the simulation. Figure 2 shows that the components of the temperature calculated across the box are constant and agree with the value dictated by eq 6. There are some small oscillations associated with the localization of the wall particles by the restraining springs. The temperature profiles calculated in the three Cartesian directions are equivalent, indicating the correct equipartition of energy among the translational degrees of freedom of the DPD particles. The elements of the pressure tensor can be calculated using the Irving-Kirkwood32,33 definition for a planar geometry

p(z) ) F(z)TI rijrij dU 1 1

〈∑

LXLY

i 0 and zero otherwise. The contribution to the pressure profile from the tethering potential associated with the wall particles is calculated in exactly the same way. Specifically, the cell is divided into slabs of width ∆z ) 0.04. The particles i and j contribute to the pressure in a particular slab if the line joining i and j, crosses, starts, or finishes in that slab. The contribution from a pair ij is distributed uniformly along the line joining the particles. Figure 3 shows the normal pressure, the difference between the normal and the tangential pressure and one off-diagonal element of p. At mechanical equilibrium, the normal pressure, pzz(z) is constant across the box reflecting the symmetry of the system. There are oscillations associated with the tethered (32) Walton, J. P. R. B.; Tildesley; D. J.; Rowlinson, J. S. Mol. Phys. 1983, 48, 1357. (33) Walton, J. P. R. B.; Tildesley, D. J.; Rowlinson, J. S. Mol. Phys. 1986, 58, 1013.

particles in the wall and they can be removed by switching from the Irving-Kirkwood method for the pressure calculation to the method of planes introduced by Todd et al.34 pN(z) - pT(z) is almost zero over the whole of the fluid region of the slab indicating that the surface tension associated with the polymer fluid interface is close to zero. The oscillations at the solid fluid interface that make it difficult to estimate σSL for this system although, we expect it to be small. The observations that pN(z) is constant across the box and the off-diagonal element, pxy(z) ) 0, indicate that the system is at mechanical equilibrium.35 We have reproduced this behavior for chains up to Nb ) 20. The structure of the fluid in the direction normal to the surface is characterized by calculating the density profile of the DPD beads

F(z) )

1

N

∑Hn(zi)

LxLy∆z i)1

(14)

Figure 4 shows the density profile for beads associated with the polymer chains for the polymers with 20 beads at coverages of 1/3 and 1/6. The solvent particles are excluded, and the only wall particles included are those (34) Todd, B. D.; Evans, D. J.; Daivis, P. J. Phys Rev. E 1995, 52, 1627. (35) Rowlinson, J. S.; Widom, B. Molecular Theory of Capillarity; Clarendon: Oxford, U.K., 1989.

4736



between brushes can be quantified in terms of I(Lz) given by

I(Lz) )

Figure 4. F(z) for the two polymer brushes with bead-length 20 at coverages of 1/3 (s) and 1/6 (‚‚‚)

z

z

z

(15)

where F′(z) is the density profile of one of the two brushes. Values of I(Lz) are given in Table 1. As expected, at fixed Lz, the overlap increases with chain length. At fixed Lz and NbI(Lz) increases with increasing coverage; e.g. Nb ) 20 and Lz ) 18.78, I(Lz) ) 0.013 for a coverage of 1/3 and 0.0019 for a coverage of 1/6. The denser the chains, the stronger the overlap of the brushes at a fixed separation. The density profiles for the solvent particles and the polymer beads as a function of the coverage are displayed in Figure 6a-d for the four different polymer chain lengths. In Figure 6a for 20 beads, we see that the solvent particles density decreases slowly from a maximum in the center of the box and remains at nonzero value approaching the wall. As one would expect for a good solvent in all the simulations, the solvent penetrates into the region of the polymer chains close to the wall. The interpenetration between the brushes vanishes as the polymer chain length decreases below 20 beads and a plateau corresponding to the pure solvent develops in Fsolvent(z). For a polymer brush, the width of the layer, w, can be defined as the distance from the wall at which the Fpolymer(z) reaches 1% of its maximum value. The widths of the layers are included in Table 1, the layers become more compact with decreasing surface coverage, and the effect is more dramatic for the longer chains. The solvent particles behave as a liquid in the center of the box between the brushes. We have calculated the in-plane radial distribution function

Figure 5. F(z) as a function of (z/Nb)1/2 for polymers of different chain lengths: Nb ) 10, 15, and 20. The curves were obtained at a coverage of 1/3 and have been displaced upward with increasing Nb for clarity. Key: (a) Nb ) 20 beads, (b) Nb ) 15 beads, and (c) Nb ) 10 beads. The lines resulting from a regression analysis are also plotted.

that form the head of the chain. Both curves are almost perfectly symmetrical which is a feature of the rapid equilibration in the DPD method. The tethering of particle 1 gives rise to four distinct peaks in F(z) for the polymer particles closest to the wall. The peaks are sharper at the higher surface coverage since the stronger in-plane correlation between the polymer chains enhances the translational ordering in the z direction. SCF theory for single brushes in a good solvent5 predicts a parabolic shape for the density profiles. We can test this prediction by plotting F(z) against (z/Nb)2 and looking for a linear correlation. Figure 5 shows F(z) against (z/Nb)2 for Nb between 10 and 20. For Nb ) 10, the brushes are well-separated and the plot is linear with a squared correlation coefficient of 0.992 between 0.16 < (z/Nb)2 < 0.63; the brush profile is parabolic away from the wall. At Nb ) 15 and 20, the plots are less linear as the brushes begin to overlap (the overlap of the brushes is quantified in Table 1). The plot for Nb ) 5 does not fit conveniently on the same scale, and the chains are so short that we might expect the profile to be dominated by wall effects. However, even in this case, F(z) against (z/Nb)2 for 0.52 < (z/Nb)2 < 0.70 is linear with a squared correlation coefficient of 0.998. There is only a small overlap between the brushes in the center of the box for the longest chains. The overlap

∫LL/2F′(z) dz/∫0L F′(z) dz

〈 g(s) )

Hn(zi)Hn(zj)δ(s - sij)〉 ∑i ∑ j*i

Hn(zi)Hn(zj)2πs ds F(zn)∆z〉 ∑i ∑ j*i

(16)

〈

where s is the in-plane separation of DPD beads. The resulting g(s) for -0.53 < z < -0.33 (i.e., for one slab in the center of the box) is typical of those observed for the isotropic DPD fluid31 but quite different from the lowdensity result for the DPD potential

a gLD(s) ) exp(-uDPD(s)/T), udpd(s) ) (r - 1)2r < 1 2 (17) The components of the mean-square displacement of the solvent particles in two different regions of the box were calculated using eqs 18 and 19. 2

∆r| (t) )

1 Ns

Ns

〈(xi(t) - xi(0))2 + (yi(t) - yi(0))2〉 ∑ i)1

∆r⊥2(t) )

1 Ns

(18)

Ns

〈(zi(t) - zi(0))2〉 ∑ i)1

(19)

These correlation functions were estimated using one time origin. They were calculated for two regions: with the solvent molecules in the solvent region -7 < z < +7 and the solvent molecules in the wall/polymer region z < -7


Langmuir, Vol. 16, No. 10, 2000 4737

Figure 6. F(z) for the grafted polymer chains and the solvent at two coverages: 1/3 coverage, solvent (s), polymer (‚‚‚); 1/6 coverage, solvent (- -), polymer (-‚‚‚-). Key: (a) Nb ) 20 beads, (b) Nb ) 15 beads, (c) Nb ) 10 beads, and (d) Nb ) 5 beads.

and z > 7. Ns is the number of solvent particles in a particular region at t ) 0. This same set of particles is followed at all later times. Parts a and b of Figure 7 show the perpendicular and parallel components of the mean square displacement for a polymer chain of length 20 at a coverage of 1/3. The diffusion in the direction perpendicular to the surface is small and liquidlike as one would expect for a confined fluid and almost independent of whether the molecule begins in the region close to the wall or in the center of the slab. The diffusion along the pore axis is significantly greater for the solvent particles in the middle of the slab than for those in the polymer/wall region. The solvent molecules are trapped in the polymer/wall region between the polymer chains reducing their displacement in the plane of the wall. In fact the diffusion coefficient in the polymer/wall region is quite isotropic. Additional information on the conformation of the polymer chains themselves can be deduced from the radii

of gyration36,37

〈S| 〉 ) 2

〈 ∑[

1

〈S⊥2〉 )

〈 [

1

N

]〉

Nb

∑ 〈(xk - xcom,i)2 + (yk - ycom,i)2〉

N i)1 Nb k ) 1 1

N

1

Nb

]〉

∑ ∑ 〈(zk - zcom,i)2〉 N i)1 N k ) 1 b

〈Stot2〉 ) 〈S|2〉 + 〈S⊥2〉

(20) (21) (22)

where the subscript “com” refers to the center of mass of a particular chain i and the subscript k refers to the beads in a particular chain. The mean squared radius of gyrations (36) Flory, P. J Statistical Mechanics of Chain Molecules; Interscience Publishers: New York, 1986. (37) Doi, M; Edwards, S. F. The theory of polymer dynamics; Clarendon Press: Oxford, U.K., 1986.

4738



Table 2. Parallel Components S|2, Perpendicular Components S⊥2, and Square Radius of Gyration Stot2 Computed during the Simulations over 100 000 Time Steps with the Standard Deviation at Two Coverages 1/3 and 1/6 for Four Different Polymer Chain Lengths coverage 1/3

coverage 1/6

polymer chain length

Stot2

S|2

S⊥2

Stot2

S|2

S⊥2

20 15 10 5

5.53 ( 0.15 3.64 ( 0.11 2.00 ( 0.06 0.77 ( 0.02

1.65 ( 0.10 1.21 ( 0.06 0.81 ( 0.05 0.37 ( 0.02

3.88 ( 0.12 2.43 ( 0.10 1.19 ( 0.05 0.40 ( 0.02

4.45 ( 0.22 2.94 ( 0.15 1.74 ( 0.09 0.74 ( 0.03

1.80 ( 0.15 1.32 ( 0.10 0.85 ( 0.07 0.40 ( 0.03

2.65 ( 0.17 1.62 ( 0.11 0.89 ( 0.07 0.34 ( 0.03

Figure 8. Mean square radius of gyration of polymer chains as a function of the number of polymer beads Nb: (a) the mean square radius of gyration parallel to the wall;(b) the mean square radius of gyration perpendicular to the wall for the two coverages 1/ (0) and 1/ (9). The inset shows how the root-mean-square 3 6 radius of gyration scales with (Nb - 1) for the two coverages 1/ (0) and 1/ (9). 3 6 Figure 7. Components of the mean square displacement as a function of time for Nb ) 20 beads at a coverage of 1/3. (a) Perpendicular component. (b) Parallel component: (s) solvent particles in the solvent region; (‚‚‚) solvent particles in the polymer/wall region.

and its components are reported in Table 2 and displayed in parts a-c of Figure 8 as a function of the number of beads in the polymer chain. The shape of the polymer is sensitive to the surface coverage. As we expected, the parallel component of the mean square radius of gyration increases as the coverage is reduced whereas the perpendicular component decreases. The inset to the figure shows how Stot scales with Nb. The results are only available over a short range of Nb, but it appears that the radius of gyration varies with the number of beads to the powers 0.63 ( 0.02 and 0.57 ( 0.02 for coverages of onethird and one-sixth, respectively. It is likely that if we extended the range of Nb, that we would find that the

exponent is coverage independent. The experimental value for the exponent of S, for a polymer in dilute solution of a good solvent, is 0.59.36,37 The values obtained in a recent DPD simulation of a polymer in solution are 0.51 ( 0.01 (weak spring model) and 0.52 ( 0.01 (strong spring model).29 The shape of the polymer can also be estimated from the tensor

SRβ )

1

Nb

∑(ri - rcom)R(ri - rcom)β

Nb i)1

(23)

which is closely related to the inertia tensor. Diagonalization of S results in three eigenvalues λ, which sum to 〈Stot2〉, and the largest of which corresponds to an eigenvector representing the long axis of the polymer. The


Langmuir, Vol. 16, No. 10, 2000 4739

asphericity38 of the polymer is defined as

A)

(λi - λj)2 ∑ i>j 3

(24)

λi)2 ∑ i)1

2(

The value of A averaged over all polymers and over time is given Table 1. If A tends to zero, the chain adopts a spherical conformation, while A tends to 1 as the chains become rodlike. For spring-bead polymers with 5 and 20 beads in a line, A takes the values 0.85 and 0.997, respectively. The values in Table 1 show that the chains stretch as the coverage increases. For Nb ) 20 at a coverage of 1/6, the asphericity parameter corresponds to a linear chain conformation with between three and four beads. The orientation of the polymer chain can be estimated from the vector joining the head of the molecule to the tail bead. (Alternatively, the orientation of a particular polymer can be associated with the principal axis obtained from the diagonalization of S. We have used both definitions and find little difference in the tilt distribution functions). Defining the chain orientation to point away from the wall, the normalized orientational distribution n(cos θ) is calculated where θ is the angle between the chain and the surface normal. n(cos θ) ) 0 for -1 e cos θ e 0 since the chains do not point back into the surface region. The chain vector can also be used to calculate the azimuthal distribution n(φ). These distributions are shown for a variety of polymer chain lengths at a coverage of 1/3 in Figure 9. As expected the azimuthal distribution is constant within the noise. The average tilt of the layer can be expressed as an order parameter with respect to the layer normal.

〈P1,n〉 )

∫-11 cos θ n(cos θ) d(cos θ)

(25)

The values of 〈P1,n〉 are given in Table 1. As the chains become longer the degree of ordering with respect to the layer normal increases at a given surface density. The ordering with respect to the surface normal decreases with decreasing coverage at fixed chain length. It also possible that the layer aligns with its own director and this can be measured using a nematic order parameter 〈P2,d〉. This is most readily achieved by calculating the Q tensor for the unit vectors ei describing the molecular axes

Q)

1

1

∑eiei - 3I N i

(26)

The largest eigenvector corresponding to the largest eigenvalue is the director d. An estimate of the tilt of the director with respect to the surface normal is given by θdirector ) cos-1(dz). The director for the layer estimated in this way is always within 10° of the surface normal. The order parameter of the polymers with respect to the director is given by

3 〈P2,d〉 ) λmax 2

(27)

where λmax is the largest eigenvector of Q. Again, this order parameter, see Table 1, increases with increasing polymer chain length at fixed coverage and decreases with decreasing coverage at fixed chain length. The orienta(38) Rudnick, J.; Gaspari, G. J. Phys. A 1986, 19, L191.

Figure 9. (a) Distribution of tilt angles with respect to the surface normal, n(cos θ). (b) Azimuthal distribution function n(φ). For different chain lengths (s) Nb ) 20 beads, (- - -) Nb ) 15 beads, (-‚-‚-) Nb ) 10 beads, (‚‚‚) Nb ) 5 beads.

tional ordering within the layer matches that of the layer with respect to the surface normal. IV. Conclusions We have demonstrated that the DPD method can be used to simulate grafted polymer brushes in a pore geometry. The temperature and pressure profiles show that the system is in thermodynamic and mechanical equilibrium. The profile of the polymer and solvent are perfectly symmetrical, indicating that the system relaxes on the time scale of the simulation. Typical simulations are of order 105 steps, which corresponds with a total real time of 2.5 ms. A total of 105 steps can be completed on a 500 MHz PC in 1 day, and the method compares favorably with the molecular dynamics method using steeply varying repulsive potentials where simulation of 2 × 106 steps may still result in poorly equilibrated profiles.39 The polymer density profiles show a parabolic shape in agreement with the SCF theory for brushes in a good solvent. The fit to a parabolic form is better for larger chain lengths, and deviations are most apparent in the overlap region between the brushes. At fixed wall separation and chain length, the overlap of the brushes increases with increasing coverage. The denser the chains, the stronger the overlap of the brushes at a fixed separation. The widths of the polymer brushes decrease significantly with decreasing surface coverage and the effect are more dramatic for the longer chains. The diffusion of the solvent beads in the direction perpendicular to the surface is small and liquidlike as one (39) Kong, Y. C.; Tildesley, D. J. Mol. Phys. 1997, 92, 7.

4740


would expect for a confined fluid and almost independent of whether the molecule begins in the region close to the wall or in the center of the slab. The diffusion along the pore axis is significantly greater for the solvent particles in the middle of the slab than for those in the polymer/ wall region. The shape of the polymer is sensitive to the surface coverage. As we expect the parallel component of the mean square radius of gyration increases as the coverage is reduced whereas the perpendicular component decreases. We have shown that the 〈Stot2〉 scales with the chain length with an exponent that is larger than that observed for the same polymer model in the bulk. Calculation of the average asphericity parameters from the inertia tensor of the polymers also shows that the chains stretch with increasing coverage. The orientational ordering of the chains is characterized by a tilt distribution, n(cos θ) that is peaked around θ ) 0, where θ is the angle between the chain axis and the surface normal. As the chains become longer the degree of ordering with respect to the layer normal increases at a given surface density. The ordering with respect to the


surface normal decreases with decreasing coverage at fixed chain length. The director of the layer has also been calculated, and the second-rank order parameter that describes the degree of ordering also increases with increasing chain length. The director is close to the layer normal. These simulations establish the DPD method as a sensible tool to study equilibrium properties of polymer brushes. It will also be possible to study the lubrication of surfaces coated with polymers by performing nonequilibrium simulations. It is possible to shear the system by moving the particles in the top and bottom walls at a fixed velocity in opposite directions and to calculate the friction from the off-diagonal elements of the pressure tensor. We are performing a study of this type using the DPD approach. Acknowledgment. P.M. would like to acknowledge the help of Unilever Research in completing these simulations. LA991396Z

Dissipative Particle Dynamics Simulations of Grafted Polymer Chains

Recommend Documents