7274
J. Phys. Chem. B 2010, 114, 7274–7285
Mesoscale Modeling of Polyelectrolyte Brushes with Salt Cyrille Ibergay,† Patrice Malfreyt,*,† and Dominic J. Tildesley‡ Laboratoire de Thermodynamique et Interactions Mole´culaires, Clermont UniVersite´, UniVersite´ Blaise Pascal, CNRS, UMR 6272, LTIM, F-63177 Aubiere, BP 10448, F-63000 Clermont-Ferrand, and UnileVer Research, Port Sunlight, Bebington, Wirral CH63 3JW, U.K. ReceiVed: December 7, 2009; ReVised Manuscript ReceiVed: April 19, 2010
We report dissipative particle dynamics (DPD) simulations of a polyelectrolyte brush under athermal solvent conditions. The electrostatic interactions are calculated using the particle-particle particle-mesh (PPPM) method with charges distributed over the particles. The polymer beads, counterions, co-ions, and solvent particles are modeled explicitly. The DPD simulations show a dependence of the brush height on the grafting density and the charge fraction that is typical of the nonlinear osmotic brush regime. We report the effect of the addition of salt on the structural properties of the brush. In the case of a polyelectrolyte brush with a high surface coverage, the simulations reproduce the transition between the nonlinear osmotic brush regime where the thickness of the brush is independent of the salt concentration and the salted regime where the brush height decreases weakly with the salt concentration. Introduction Polymers grafted onto surfaces have applications in many important technical problems, wetting, adhesion, colloid stabilization, and biocompatibility. At high surface coverages in good solvent conditions, the polymer chains are strongly stretched in the direction perpendicular to the grafted surface, giving rise to a brush structure. In the case of a neutral polymer brush, this structure results from a balance between the entropic elastic energy of the chain, the monomer-monomer interactions, and the affinity of the polymer for the solvent. In the case of a charged polymer, a polyelectrolyte brush, the structure is influenced by additional long-range electrostatic interactions that depend on the fraction of monomers carrying an ionic charge and the concentration of salt in solution. The brush thickness is then determined by a balance between the osmotic pressure of the confined counterions and entropic elasticity of the chains opposing the extension. The dependence of some structural properties of the polyelectrolyte brushes such as the brush height with respect to the grafting density (Fa), the degree of ionicity (f), and the added salt concentration (cs) has been previously studied by theory,1-6 simulations,5-12 and experiment.13,14 Pincus2 has shown that the physical behavior of polyelectrolyte brushes can be predicted by scaling theories. In the osmotic brush regime, all of the counterions are trapped in the brush layer if the Gouy-Chapman length λGC is small compared to the brush height h. In this regime, scaling arguments suggest that the brush height is independent of the grafting density. However, a number of simulations5,9,11 and experiments6,14 have shown that the brush thickness increases with the surface coverage. To understand these results, a nonlinear osmotic brush regime5,6,11 has been suggested. In this regime, the brush is in the strong-charging and strong-stretching limits, and the brush height is found to scale with the grafting density as ∼(f + d2Fa)/(1 + f), where d is a measure of the monomer and counterion diameters. When * To whom correspondence should be addressed. E-mail: Patrice.
[email protected]. † Clermont Universite´, Universite´ Blaise Pascal. ‡ Unilever Research.
salt is added, the system follows the brush osmotic regime, and the brush collapses with increasing salt concentration with a weak power law of h ≈ cs-1/3.2 Alternative theories15 predict h ≈ cs-2/3. Previous simulations of polyelectrolyte brushes have been performed using both the Monte Carlo (MC) and molecular dynamics (MD) methods using coarse-grained potentials. Monte Carlo off-lattice simulations of polyelectrolytes7 were performed without the explicit inclusion of counterions or the excluded volume interactions. More recently, MC simulations have been developed to study the property of short charged polymers grafted to surfaces.12 In this case, the electrostatic interactions were calculated using the Ewald method, and the counterions were included explicitly. In the MD simulations of polyelectrolyte brushes,5,6,8-11 the polymer chains were constructed using a bead-spring model. The solvent was treated in the primitive model, using a Bjerrum length for water of 7.1 Å. Good solvent conditions were obtained by using a purely repulsive short-range Lennard-Jones (LJ) interaction between monomers and monomers and counterions. The use of the steep LJ potential requires a relatively small time step for integrating the equations of motion and an additional white noise force, familiar in Brownian dynamics, which is added to improve the thermal equilibration. The solvent, beyond the counterions, was not treated explicitly in these simulations. A genuinely mesoscale approach to modeling the brush is to use the dissipative particle dynamics (DPD) technique proposed by Hoogerbrugge and Koelman.16,17 DPD is a stochastic dynamics method in which the individual particles (beads) represent a set of atomistic degrees of freedom (that is, a region of the solvent or the polymer). As a consequence of this coarsegraining, the conservative interactions are soft. The dissipative and Brownian forces are short-ranged and pairwise-additive, so that the hydrodynamics of the fluid satisfies the Navier-Stokes equation. There are three advantages to the DPD method. First, the soft potential allows for a time step that is up to an order of magnitude larger than those typically used in MD simulations. Second, the method, unlike the Brownian dynamics technique,
10.1021/jp9115832 2010 American Chemical Society Published on Web 05/10/2010
Mesoscale Modeling of Polyelectrolyte Brushes with Salt includes a correct description of hydrodynamic interactions. Third, the solvent can be explicitly included by modeling regions of the solvent as individual particles. It is also possible to modify the standard DPD algorithm to prevent artificial bond crossing of the soft polymers18-20 and to include distributed charge distributions to model polyelectrolytes.21-23 In the past, we have successfully applied DPD to study the equilibrium properties18,24,25 and rheological properties20,26,27 of neutral polymer brushes. The DPD methodology was modified to allow the shearing26 of polymer brushes. The friction coefficient was calculated as a function of the solvent quality, and a correlation between the interpenetration coefficient and the frictional forces was established. DPD simulations in the grand canonical ensemble18,20,25,27 were able to reproduce the main structural and rheological properties of the neutral brushes under compression. The modeling of polyelectrolyte brushes requires the calculation of the long-range electrostatic interactions at a mesoscopic level. Groot21 has proposed their incorporation using an adapted version of the particle-particle particle-mesh (PPPM) method. Gonzales-Melchor22 has proposed a method combining the Ewald28 technique with a charge distribution on the DPD particles. We have tested these two approaches23 and have shown that the structural and mechanical properties of the polyelectrolyte brushes using both techniques are essentially equivalent. However, the PPPM method is faster than the Ewald summation; PPPM scales as O(N log N), whereas the Ewald method scales as O(N3/2), where N is the number of charged particles. Recently, the DPD method associated with the PPPM technique has been successfully applied to investigate the overcharging behavior of a cylindrical polyelectrolyte brush in the presence of multivalent counterions29 and the interactions of the complex between a cylindrical polyelectrolyte brush and a linear polyelectrolyte with the opposite charges.30 The method has also been applied to the study of the interaction31 between two polyelectrolyte brushes in aqueous media. The brushes were formed from diblock copolymers composed of a short hydrophobic anchoring block and a long hydrophilic block. The thrust of our research is to study the interaction between polyelectrolyte brushes under shear, but in this paper, we focus on the properties of a single polyelectrolyte brush in a salt solution as predicted by the DPD method. We use a standard set of DPD parameters21 to model grafted homopolyelectrolyte chains. We will address the following questions: how does the explicit inclusion of the solvent by particles affect the structure of the polyelectrolyte brush? Does the DPD model, with its soft repulsive potential, reproduce the scaling behavior of the nonlinear osmotic regime? In this paper, we will study the properties of a single polyelectrolyte brush as a function of the surface coverage and charge fraction. In particular, we will study the dependence of the brush height on Fa and f and compare with the scaling laws of the nonlinear osmotic regime. We complete the study by the modeling of polyelectrolyte brushes in the presence of excess salt. We will investigate the impact of the salt on the structural properties of the brush and compare our mesoscale simulations with previous MD simulations5,6,8,9,11 and theoretical predictions.1-6 Dissipative Particle Dynamics (DPD) Model Standard DPD Forces. In the DPD approach, solvent particles are coarse-grained into soft beads that interact with the pairwise additive force fi defined as the sum of three contributions
J. Phys. Chem. B, Vol. 114, No. 21, 2010 7275
∑ (fijC + fijR + fijD)
fi )
(1)
j*i
where fCij , fRij , and fDij are the conservative, random, and dissipative forces, respectively. The conservative repulsive force fCij is derived from a soft interaction potential and is expressed as follows
fijC )
{
aijωC(rij)rˆij (rij < rc) (rij g rc) 0
(2)
where aij is the maximum repulsion parameter between particles i and j, rij is the relative displacement of particles i and j, and rˆij is the corresponding unit vector. rc is the cutoff radius. The weight function ωC(rij) is equal to 1 - rij/rc for rij e rc and vanishes for rij g rc. The dissipative and random forces are given by
fijD ) -γωD(rij)(rˆij.vij)rˆij fijR ) σωR(rij)θij
1 rˆij √δt
(3)
(4)
where δt is the time step. vij ) vi - vj is the relative velocity, σ is the amplitude of the noise, and θij is a random Gaussian number with zero mean and unit variance. γ and σ are the dissipation strength and noise strength, respectively. The terms ωD(rij) and ωR(rij) are dimensionless weighting functions. Espanol and Warren32 have shown that the system will sample the canonical ensemble and obey the fluctuation-dissipation theorem if the following conditions are satisfied
γ)
σ2 2kBT
and ωD(rij) ) (ωR(rij))2
(5)
where kB is Boltzmann’s constant and T is the temperature. The weighting function ωR(rij) is chosen to be similar to ωC(rij). The equations of motions are integrated using a modified version of the velocity-Verlet algorithm.21 The force is updated once per iteration, and because the force depends on the velocities, the velocity in the next time step has to be estimated by a predictor algorithm. The velocity is then corrected in the last step. The reduced time step δt was taken to be 0.06 for all of the simulations reported in the paper. When modeling polymers, the integrity of the chain is ensured by an additional spring force between neighboring beads given by
fijS ) -ks(rij - r0)rˆij
(6)
where the equilibrium bond distance r0 is zero and the spring constant ks ) 4.0. This pairwise force21 is then added to the sum of the DPD conservative force of eq 1. Electrostatic Interactions (PPPM Method). We use an adapted version of the PPPM method for the calculation of the electrostatic forces within the DPD methodology. The electrostatic interactions were originally incorporated into the DPD model by Groot.21 The first step consists of finding a model for the density distribution adapted to a charged bead. The use of
7276
J. Phys. Chem. B, Vol. 114, No. 21, 2010
Ibergay et al.
a soft potential in DPD allows the overlap between DPD beads. When fully charged DPD beads are modeled, this can lead to the formation of artificial ion pairs and cause the divergence of the electrostatic potential at r ) 0. To avoid this problem, Groot chose to spread out the charges using the following distribution
f(r) )
3 (1 - r/re) πr3e
for r < re
(7)
where re is the electrostatic smearing radius; f(r) ) 0 when r is greater than re. The interaction potential between two smearedout charges is given in eq A-1. The potential and the corresponding force are represented in Figure A-1 of the Appendix A. The value of the smearing radius re is justified in Appendix A. The electrostatic field is then solved on a lattice according to the method of Beckers et al.33 The charges are assigned to the lattice nodes within the cell, and the long-range part of the interaction potential is calculated by solving the Poisson equation on the grid. Details of the charge assignment can be found in Groot’s21 original paper. Whereas in the original PPPM method the electrostatic field was solved by using a Fourier transform, the method developed by Groot used a real-space successive overdamped relaxations method. This makes the Groot’s method close to the multigrid method of Sagui and Darden.34 The electrostatic force fEi on a charged bead i is calculated from
fiE(ri) ) -qi
∑ fj(ri)∇Ψ(rj)
(8)
j
where ri is the position of the charged bead i and qi is the number of unit charges on bead i. fj(ri) is defined as
fj(ri) )
1 - |rj - ri |/re
∑ (1 - |rn - ri|/re)
(9)
n
where rj is the position of the node j, and the sum over n runs over all nodes with a distance re from ri. This function means that a charge proportional to f(r) in eq 7 is assigned to each node i and normalized such that the sum of all of the charges within a distance re equals the charge on the bead i. Ψ(rj) is the local electrostatic field at lattice node j. The field Ψ(rj) is calculated from the Poisson equation expressed in reduced DPD units as
∇* · (p(r)∇*Ψ) ) -ΓF*
(10)
where F* is the concentration of cations minus that of anions per r3c , ∇* is the gradient in reduced units, and p(r) is the local polarizability relative to pure water. The total momentum of the simulation cell is conserved by removing for each charge a small residual force. This residual force is on the order of 5 × 10-5 in reduced units and is expressed as ΣNi qifEi /N. The total force fi of eq 1 is then modified by adding the fEi contribution and the residual force. In the simulation, the particle mass, temperature, and interaction range are chosen as units of mass, energy, and length; hence, m ) kBT ) rc ) 1.0. The unit of time τ becomes then rc(m/ kBT)1/2. The real length rc can be estimated from the volume of DPD bead, and then, rc ) (F*NmVm/NA)1/3, where F* is the
reduced density of DPD particles, Vm ) 18 cm3 mol-1, and NA is Avogadro’s number. Groot21,35 used Nm ) 3 and a reduced density of 3. The coupling constant Γ is given by e2/(kBTε0εrrc), where e is the electronic charge, ε0 ) 8.85418782 × 10-12 C2 J-1 m-1 is the permeability of free space, and εr ) 78.3 the relative permittivity of water at 298 K. Using the previous values, rc ) 6.46 Å and Γ ) 13.87. As already discussed by Groot and Rabone,35 the time scale is fixed by matching the diffusion constant of water. For a repulsion parameter a ) 25, we find that the natural unit of time τ is 160 ps. The Bjerrum length, defined by λB ) e2/(4πkBTε0εr), characterizes the length scale at which two unit charges interact with thermal energy kBT in a medium of dielectric constant ε. With the parameters used in this work, λB is then equal to 1.11 in reduced DPD units. The thermal and mechanical equilibria of polyelectrolyte brushes with salt are demonstrated through the calculation of the profiles of the kinetic and configurational parts of the pressure tensor calculated using the method of planes36 (see Figure B-1 of Appendix B). Note that the electrostatic force fEi is not pairwise-additive and cannot be used within the Irving and Kirkwood formalism37 to calculate the pressure tensor. Simulation Cell The system consists of two planar solid surfaces composed of three layers of 324 DPD particles. The two surfaces are positioned at the top and the bottom of the simulation cell. The presence of the second surface allows one to obtain finite values of concentration due to the finite length along z. One of the two surfaces is coated with Np polymer chains, which are randomly grafted by a harmonic force acting between the end particles of the chains and the particles of the first layer of the wall. Each chain contains Nb ) 20 or 30 polymer beads. The surface coverage can defined by Fa ) Np/(LxLy). Different charge fractions f are studied, 1.0 (completely negatively charged), 0.5 (half fully charged), and 0 (neutral). To preserve the neutrality, there are N+ ) f × Np × Nb counterions. The number of solvent particles Ns is adjusted so that the overall reduced density between the two walls is close to 3.0. The aij parameters are set to 25.0 for all interactions. The polymer brush is then modeled in athermal solvent conditions. The cell dimensions are Lx ) 16.7rc and Ly ) 6.4rc. The Lz dimension of the primary cell is Lz ) 53 in reduced units, corresponding to a separation distance D between the two surfaces of 50. To respect the supercell approximation, the simulation cell was elongated along the z direction by placing an empty space of at least twice the space of the fluid-occupied region. The periodic boundary conditions were then applied in the three directions. The resulting Lz′ was fixed to 152. The simulations consisted of an equilibration period of 100 000 steps, followed by an acquisition period of 300 000 steps. Results and Discussions Effect of the Surface Coverage. Figure 1a shows density profiles of solvent, counterions, and monomers as a function of the distance from the surface for a grafting density of 1.01 and for completely charged chains (f ) 1.0). We observe that the solvent particle density decreases slowly from a maximum in the bulk-like region of the box and remains at a nonzero value approaching the wall. As one would expect for an athermal solvent, the solvent penetrates into the region of the polymer chains close to the wall. The monomer density is greater than that of the solvent at a z value corresponding to half of the brush height. This figure suggests that the explicit presence of solvent
Mesoscale Modeling of Polyelectrolyte Brushes with Salt
J. Phys. Chem. B, Vol. 114, No. 21, 2010 7277 TABLE 1: Brush Height 〈h〉, Thickness of the Counterion Layer 〈zci〉, Mean Distance between Two Connected Beads 〈b〉, and Degree of Condensed Counterions 〈fci〉 as a Function of the Surface Coverage Defined in Terms of s and Ga ) Np/LxLya
Figure 1. Monomer, counterion, and solvent density profiles for a polyelectrolyte brush. The polymer chain length is fixed to 20 beads, and the surface coverage is Fa ) 1.01. The solvent, monomer, and couterion profiles are shown on the same scale in (a). Panel (b) focuses on the density profiles of the monomers and counterions.
molecules inside of the brush and outside of the brush may lead to differences in the shape of monomer density profiles with simulation methods that do no treat the solvent molecules explicitly. Additionally, Figure 1b shows that the profiles of the polymer beads and counterions are coincident, indicating that the counterions are almost completely confined in the brush layer. Note that the precise degree of confinement depends on both the degree of charging and surface coverage. We begin by investigating the impact of the surface coverage on the structural properties of the polyelectrolyte brushes. The grafting density is changed from a small value of Fa ) 0.06 with only six polymer chains grafted to the surface to a relatively high value of Fa ) 1.01 where one-third of the wall particles of the outermost layer are connected to a polymer bead. The number of grafted polymer chains and solvent particles is given for each surface coverage in Table 1. Figure 2a shows the monomer density profiles of a fully charged brush at various grafting densities. The polymer chain length is fixed at 20 beads. This figure shows that the profiles exhibit the same features at each grafting density, a strong ordering close to the grafted surface due to the wall effects and an increase of the elongation of the brush along z with the grafting density. Indeed, the balance between the tethering of the first bead of the polymer chain to the surface and the repulsive interactions between the polymer beads and the wall particles gives rise to three distinct peaks in F(z) for the polymer beads close to the wall. The peaks are sharper at the higher surface coverage since the stronger in-
s
Fa
Np
〈b〉
〈h〉
〈zci〉
1/3 1/4 1/5 1/6 1/7 1/8 1/20 1/30 1/40 1/50
1.01 0.76 0.61 0.51 0.43 0.38 0.15 0.10 0.07 0.06
108 81 65 54 46 41 16 11 8 6
f ) 1.0, Nb ) 30 9543 3240 11163 2430 12123 1950 12783 1620 13263 1380 13563 1230 15063 480 15363 330 15543 240 15663 180
1.15 1.11 1.09 1.09 1.08 1.07 1.06 1.06 1.06 1.06
28.1 25.5 24.0 23.1 22.3 21.9 19.5 19.4 19.0 18.6
28.3 26.3 24.7 23.9 23.0 22.4 20.2 19.5 19.5 19.3
1/3 1/4 1/5 1/6 1/7 1/8 1/20 1/30 1/40 1/50
1.01 0.76 0.61 0.51 0.43 0.38 0.15 0.10 0.07 0.06
108 81 65 54 46 41 16 11 8 6
f ) 1.0, Nb ) 20 11703 2160 12783 1620 13423 1300 13863 1080 14183 920 14383 820 15383 320 15583 220 15703 160 15783 120
1.14 1.11 1.09 1.08 1.08 1.07 1.06 1.06 1.07 1.07
19.3 17.7 16.7 16.0 15.6 15.5 14.0 13.9 13.5 13.2
20.2 18.2 16.8 16.4 16.3 15.8 14.5 14.5 14.0 14.0
1/3 1/4 1/5 1/7 1/10 1/30 1/50
1.01 0.76 0.61 0.43 0.30 0.10 0.06
108 81 65 46 32 11 6
f ) 0.5, Nb ) 20 12783 1080 13593 810 14073 650 14643 460 15063 410 15693 110 15843 60
1.06 1.03 1.02 1.01 1.00 0.99 0.98
15.5 14.2 13.4 12.4 11.9 10.4 10.3
15.8 14.5 13.7 12.8 12.3 10.9 10.8
1/3 1/4 1/5 1/7 1/10 1/30 1/50
1.01 0.76 0.61 0.43 0.30 0.10 0.06
108 81 65 46 41 32 6
f ) 0.0, Nb ) 20 13863 14403 14723 15103 15383 15803 15903
0.98 0.96 0.96 0.95 0.95 0.94 0.94
11.3 10.3 9.5 8.9 8.2 6.4 6.2
Ns
N+
a s represents the percentage of wall particles of the last layer connected to the first beads of the polymer chains. Two different polymer chain lengths Nb and two charge fractions f are used. Np, Ns, and N+ ) f × Np × Nb are the numbers of polymer chains, solvent particles, and counterions, respectively.
plane correlation between the polymer chains enhances the translational ordering along z. The number and the amplitude of these oscillations26 depend on the repulsive interaction between the polymer beads and the surface particles and on the spring constant controlling the stiffness of the chain. In Figure 2b, we plot the same profiles normalized by Nb ) ∫0∞ F(z) dz. The figure shows that the polymer chain stretches along the direction normal to the surface with increasing grafting densities. At low surface coverages, the polymer chain increasingly samples the region close to the surface, and the structure of the brush is typical of that of the mushroom regime.38,39 To measure the extension of the brush, we calculate the average height of the brush 〈h〉 and of the counterion layer 〈zci〉 from the first moment of the density profiles of monomers and counterions as
7278
J. Phys. Chem. B, Vol. 114, No. 21, 2010
∫-LL /2/2 zFm(z) dz 〈h〉 ) 2 L /2 ∫-L /2 Fm(z) dz z
z
z
z
Ibergay et al.
∫-LL /2/2 zFci(z) dz 〈zci〉 ) 2 L /2 ∫-L /2 Fci(z) dz z
z
z
z
(11) where Fm(z) and Fci(z) are the density profiles of the monolayer and the counterions, respectively. The factor of 2 takes into account for the fact that the brush height is twice the first moment when the monomer density profile is uniform inside of the brush. We use the definition suggested in previous numerical simulations.9,18,20,26 The different brush heights are reported in Table 1 with the average distance 〈b〉 between neighboring beads within the polymer chain. This table shows that the brush height is similar to the height of the counterion layer, indicating that the counterions are mainly confined to the brush layer regardless of the grafting density and the charge fraction. As observed in Figure 1, the density profiles of the monomers and counterions are very similar. The Gouy-Chapman length λGC, defined as 1/(2πλBNbfFa), is a typical length which defines the height at which counterions are effectively bound to a surface of charge density of efNbFa.2 In strongly charged
brushes, the height of the counterion layer is equal to h + 3λGC/ 2. For the strongest charged polyelectrolyte brush (f ) 1.0) and the highest grafting densities (Fa ) 1.01), the Gouy-Chapman length is very small and equal to 0.0047. The maximum value of λGC obtained with the set of charge fraction and grafting densities investigated here is 0.24. It means that λGC is quite small compared to the values of the brush height for our different systems. This confirms that all of the counterions are trapped in the brush layer and implies a local electroneutrality23 in the direction normal to the surface. Figure 2c displays the Fm(z) monomer density profiles for two chain lengths (Nb ) 20 and 30) as a function of (z/Nb) for three surface coverages. The shape of the profile is independent of the chain length. We then use a Nb ) 20 polymer for the next simulations. Panel (d) of Figure 2 shows the density profiles as a function of (z/Nb)2. The fact that the plots are linear means that the density profiles of the polyelectrolyte chains can be accurately represented by parabolic functions for the different surface coverages investigated here. The deviations observed close to the grafting surface result from the strong ordering and the packing effects close to the surface. These deviations are then more important at high surface coverages and vanish at small grafting densities. Recent density
Figure 2. Monomer density profiles as a function (a) the grafting density, (b) normalized in such a way that the integral is equal to Nb. The monomer density profiles of (c) and (d) are expressed as a function of (z/Nb) and (z/Nb)2, respectively. The solid line in (d) results from a regression analysis.
Mesoscale Modeling of Polyelectrolyte Brushes with Salt
J. Phys. Chem. B, Vol. 114, No. 21, 2010 7279 established that the brush height increased with the surface coverage. This led the authors to define a new regime called the nonlinear osmotic brush regime.5,6,41 This regime combines the high stretching of the chain elasticity with the nonlinear entropic effects of the counterions due to the finite volume of the grafted polymer chains.6 The resulting brush height is dependent on both the charge fraction and the surface coverage
f + d2Fa h = Nbb 1+f
Figure 3. Brush height h/(Nb〈b〉) as a function of the grafting density Fa for polyelectrolyte chains of Nb beads with charge fraction f. The solid lines show a fit from eq 12 for charged polymers and from F1/3 a for neutral polymers.
functional theory calculations40 also established a parabolic form for the polyelectrolyte brushes. Additionally, recent neutron reflectivity measurements14 on a polyelectrolyte brush in the salted brush regime confirm a parabolic monomer density profile. It is also then interesting to note that the monomer profiles calculated from recent MD,5,8,9 where the solvent particles are not treated explicitly, exhibit a step-like form. The inclusion of the solvent particles changes the shape of the polymer density profiles. This observation is supported by our previous DPD simulations on wet and dry neutral brushes;26 the density profile of the dry brushes (without solvent) are steplike, in contrast to the parabolic profiles associated with the wet brushes in good or athermal solvent conditions. In good solvent conditions, solvent particles penetrate significantly in the polymer-rich region and create excluded volume for polymers. The free ends of the polymer chains are distributed with nonzero density within the brush; the monomer profile is then parabolic. The wet brush exhibits a smaller brush thickness compared to that of the dry brush. This is due to the volume exclusion of the explicit solvent particles and the resulting larger normal pressure. The origin of the extension of the brush thickness in the dry brush is entropic because the interactions between the polymer beads are solely repulsive. When the wet brush is modeled in poor solvent conditions,26 no solvent particles enter the region between the polymers. As a result, the brush height decreases, and the monomer density profile has the form of a step function rather closer to that of the dry brush.FromtheseDPDsimulationsandpreviousMDsimulations5,8,9 with an implicit modeling of the solvent molecules, we can conclude that the origin of the step-like form of the monomer density profiles can be explained by the absence of explicit solvent particles within the polymer brush. In Figure 3, we show the ratio between the average brush height 〈zm〉 and the contour length defined by Nb〈b〉 for polyelectrolyte brushes with two degrees of charging (f ) 1 and 0.5) and two polymer chain lengths (Nb ) 20 and 30). For comparison, we add the ratio calculated for a neutral brush. Figure 3 shows that the brush height of the polyelectrolyte brush increases with the surface coverage, whereas the scaling theory of Pincus2 does not predict any dependence of the brush height on the grafting density. However, recent simulations5,6,9 have
(12)
where d2 is a finite effective diameter taking into account the monomer and counterion diameters. Figure 3 shows that the fully charged chains are stretched from 60 to 80% of their contour length when the grafting density is increased from 0.06 to 1.01. This elongation is reduced when the charged fraction is halved. The extension of the chains is much less pronounced for neutral polymer chains. The ratio is about 30% for low surface coverages. We need a very high surface coverage of neutral chains to recover the same extension as a polyelectrolyte chain at low surface coverages. These curves show that the charged chains are highly stretched as compared to the uncharged chains and confirms the use of the strong-stretching limit for the polyelectrolyte brushes investigated here. In Figure 3, we show that the dependence of the brush height with respect to the surface coverage can be reproduced using eq 12 for polyelectrolyte brushes with different chain lengths and fractions of charge. The parameter d is difficult to associate with a particular DPD parameter and is then treated here as a free fitting parameter. Our DPD simulations of the polyelectrolyte brushes agree very well with previous simulation results5,6,9 and with recent experiments14 and confirm an increase of the brush thickness with the grafting density. For the neutral polymer brush, the dependence of the brush height with respect to the surface coverage is found to follow the theoretical prediction38,42 h/(Nbb) = F1/3 a for relatively high surface coverages. This has been confirmed by recent computer simulation studies.18,43 Charge Fraction. In this section, we present the impact of the charge fraction of the polymer chains on the brush height. The charge fraction varies from 0.0 to 1.0 at two grafting densities (Fa ) 1.01 and 0.15). The composition of each system is given in Table 2. Figure 4a presents the monomers density profiles for Fa ) 1.01 as we move from a neutral polymer brush to a full charged polyelectrolyte brush. The density profiles extend further in the direction normal to the surface with increasing charge fraction. For examples, the brush layer is increased by a factor of 1.7 from a neutral brush to a full charged polyelectrolyte brush at a high grafting density; this factor is about 1.9 for a low surface coverage. The heights of the brush and counterion layers are given in Table 2. The results are shown in Figure 4b. The relationship between the brush height and the fraction of charge is not linear and is accurately described by eq 12. We see that the impact of increasing the charge fraction from f ) 0.0 to 1.0 on the brush height is effectively the same in magnitude as increasing the grafting density from Fa ) 0.06 to 1.01. The DPD simulations confirm the fact that the polyelectrolyte brushes follow the nonlinear brush regime for the dependence of the brush thickness with the charge fraction. The Gouy-Chapman length changes from 0.05 to 0.24 in the range of fractions of charge studied here. This means that the brush layer matches the counterion layer (see Table 2). The density profile of the counterion (not shown here) accurately
7280
J. Phys. Chem. B, Vol. 114, No. 21, 2010
Ibergay et al.
TABLE 2: Brush Height 〈h〉, Thickness of the Counterion Layer 〈zci〉, Mean Distance between Two Connected Beads 〈b〉, and Degree of Condensed Counterions 〈fci〉 as a Function of the Charge Fraction f for Two Different Surface Coveragesa 〈h〉
〈zci〉
〈fci〉
13863 13431 12999 12783 12567 12135 11700
s ) 1/3, Fa ) 1.01 0 0.97 11.3 432 1.00 13.0 864 1.04 14.7 1080 1.05 15.5 1296 1.07 16.2 1728 1.11 17.7 2160 1.15 19.3
0 13.1 14.7 15.6 16.3 17.8 19.4
0.851 0.952 0.964 0.972 0.983 0.995
15703 15639 15575 15543 15511 15447 15383
s ) 1/20, Fa ) 0.15 0 0.98 8.1 64 1.02 9.6 128 1.04 11.0 160 1.05 11.5 192 1.06 12.3 256 1.09 13.3 320 1.15 15.1
0 9.6 11.0 11.4 12.3 12.8 15.1
0.235 0.384 0.433 0.487 0.563 0.604
f
Ns
0.0 0.2 0.4 0.5 0.6 0.8 1.0 0.0 0.2 0.4 0.5 0.6 0.8 1.0
N+
〈b〉
a Np, Ns, and N+ ) f × Np × Nb are the numbers of polymer chains, solvent particles, and counterions, respectively.
matches that of the monomers and implies a local electroneutrality at all of the charge fractions. To demonstrate that the monomers are trapped in the brush layer, we calculate the polyion-counterion pair distribution function p(r), where r is the separation distance between the counterion and the closest polyelectrolyte bead. The distributions are normalized so that 2π∫∞0 rp(r) dr ) 1 and are shown in Figure 4c for three fractions of charge at two different grafting densities. These distributions clearly show that the strength of correlation between polyion and positive ions is growing with both charge fraction and grafting density. At high grafting density, the distribution functions show a pronounced peak even at a small degree of charging (f ) 0.2). At a low surface coverage, the peaks are much less pronounced and extend over larger distances. The integration of the p(r) distribution function from 0.0 to λB gives an estimate of the degree of counterion condensation within the brush layer, where λB is the Bjerrum length. With the parameters used in this work, the Bjerrum length is equal to 1.11 in reduced DPD units. The degree of counterion condensation is
∫0λ rp(r) dr ∫0∞ rp(r) dr B
〈fci〉 )
(13)
and is reported in Table 2. For a high grafting density, the degree of counterion condensation increases from 85 to 99% as the charge fraction increases from 0.2 to 1.0. In the case of a moderate surface coverage, the condensation of counterion increases from 23 to 60%. Although the degree of counterion condensation increases with grafting densities, the effect is more pronounced at lower grafting densities. The number of condensed counterions is a measure of the close association between a polymer bead and a counterion, whereas the extent of counterion confinement indicates the number of counterions trapped in the brush. The confinement changes from 96.8 to 99.9% for the range of charge fraction investigated here. All of the counterions take part in the screening of the electrostatic
Figure 4. (a) Monomer density profiles of polyelectrolyte brushes with different charge fractions f. The polymer chain length is fixed to Nb ) 20 beads; (b) brush height h/(Nb〈b〉) as a function of the charge fraction for two surface coverages. The solid lines result from the fit of the brush height with respect to f using eq 12; (c) polyion-counterion distribution function calculated for different charge fractions and grafting densities.
interactions. As a consequence, the electroneutrality is satisfied locally within the brush even at small grafting densities. Salt Concentration. In this section, we will investigate the effect of adding salt on the structure of polyelectrolyte brushes. The conservative interaction parameters for the solvent-salt, polymer-salt, and salt-salt interactions are identical (aij ) 25) to those for all the other interactions. cs is the concentration of added salt in the simulation cell of volume LxLyD, where D is the separation of the two surfaces. Figure 5a shows the profiles of the solvent particles, counterions, co-ions, and polymer beads along the surface normal for a reduced grafting density of 1.01 and cs ) 0.239. This salt concentration corresponds to the
Mesoscale Modeling of Polyelectrolyte Brushes with Salt
J. Phys. Chem. B, Vol. 114, No. 21, 2010 7281 TABLE 3: Concentrations of Free Ions Inside of the Brush csibrush and of Free Ions in the Bulk csibulk as a Function of the Salt Concentration csa csibulk
〈h〉
〈b〉
〈fci〉
2166 2224 2351 2542 2797 3116 3434 4160 5160
s ) 1/3, Fa ) 1.01 6 0.99 0.04 64 1.00 0.07 191 1.02 0.15 382 1.03 0.25 637 1.06 0.39 956 1.10 0.56 1274 1.13 0.73 2000 1.22 1.12 3000 1.35 1.62
19.3 19.2 19.0 18.9 18.6 18.3 18.1 17.7 17.2
1.15 1.14 1.14 1.13 1.13 1.12 1.12 1.10 1.10
0.987 0.989 0.997 0.998 0.998 0.999 0.999 1.000 1.000
13713 12957 12353 11599 10843 9863 7863
1155 1533 1835 2212 2590 3080 4080
s ) 1/6, Fa ) 0.51 75 0.60 0.06 453 0.66 0.25 755 0.73 0.39 1132 0.79 0.56 1510 0.86 0.73 2000 1.01 0.89 3000 1.14 1.06
15.9 15.3 14.5 14.1 13.9 13.4 12.9
1.06 1.05 1.05 1.04 1.03 1.02 1.01
0.428 0.632 0.755 0.824 0.892 0.952 0.984
14903 14303 13303 11303 7303
560 860 1360 2360 4360
s ) 1/18, Fa ) 0.17 200 0.27 0.10 500 0.34 0.23 1000 0.48 0.43 2000 0.75 0.83 4000 0.90 1.10
12.5 11.7 10.6 9.6 8.6
1.03 1.02 1.00 0.99 0.98
0.054 0.133 0.318 0.605 0.906
cs
Ns
N+
0.001 0.012 0.036 0.072 0.119 0.179 0.239 0.374 0.562
11691 11575 11321 10939 10429 9791 9155 7703 5703
0.014 0.085 0.141 0.212 0.283 0.374 0.562 0.037 0.094 0.187 0.374 0.749
N-
csibrush
a The heights of the brush and of the counterion layers are also given along with the degree of counterion condensation as a function of cs. Ns and N- are the numbers of solvent particles and added negative anions. N+ ) N- + f × Np × Nb is the number of positive ions in the simulation cell. The different properties are given for three surface coverages with Nb ) 20 and f ) 1.0.
Figure 5. (a) Monomer, counterion, and co-ion density profiles for a polyelectrolyte brush. The polymer chain length is fixed to 20 beads, the surface coverage is Fa ) 1.01, and the salt concentration cs is 0.239. (b) Local net charge as a function of the distance from the grafting surface for different salt concentrations. The profiles for cs ) 0.24 and 0.56 are shifted by 0.2 and 0.4, respectively, in the y direction. (c) Monomer density profiles as a function of z for different salt concentrations.
addition of 1274 mobile positive ions and 1274 mobile negative anions in the total volume with respect to the unsalted system (see Table 3). For this salt concentration, the co-ion density profile is homogeneous in the bulk solution and matches the counterion profile as expected for local electroneutrality in this region. At the top edge of the brush, the density profile of the co-ions starts to decrease and vanishes within the brush. No co-ion is present in the region close to the grafted surface. Local electroneutrality in the brush region implies that co-ions at the end of the polyelectrolyte chains are compensated for by the counterions. As a result, the counterion profile differs slightly
from that of the polymer beads at the brush rim. Figure 5b shows the profiles of the sum of the charged species (co-ions, polyions, and counterions) as a function of z. The electroneutrality is satisfied at each z except close to the grafted surface due to the layering of the grafted polymers. We also observe the formation of a dipole at the rim of the brush. The strength of the dipole decreases with increasing added salt concentration. Figure 5c shows the density profiles of polymer beads calculated at different salt concentrations. As the salt concentration is increased, the chains shrink in the direction normal to the surface with a swelling of the brush in the region close to grafting; the monomer density profiles maintain their parabolic shape and fall to zero more abruptly that those calculated in unsalted polymer brushes. This is in agreement with what has been observed from neutron reflectivity measurements on a monolayer brush with salt.14 Figure 6a shows the brush height h of the polyelectrolyte brush as a function of cs for three surface coverages. The brush thickness h is normalized by the brush height h0 of the corresponding unsalted brush. The curve scales as suggested by experimental13,14 and other simulation11 studies. This curve underlines two regimes for the highest grafting density studied here, the nonlinear osmotic brush regime2 where the brush height is independent of the added salt concentration and the salted regime characterized2,3,44 by a weak dependence of the brush height with respect to cs. The transition from the osmotic brush to the salted brush regime occurs when the bulk salt concentration equals the free mobile counterion concentration within the brush. Our simulations show the brush height as cs to the power of -0.15 rather than to the power of -1/3 expected from scaling theory.2,3,44 The results support the work of Kumar and Seidel,11 who observed the same power and
7282
J. Phys. Chem. B, Vol. 114, No. 21, 2010
Ibergay et al. free ions inside of the salted and unsalted brush, respectively. brush csi,0 is calculated by counting the number of free positive and negative anions within the layer defined along the z direction by the brush thickness. Interestingly, this curve exhibits the theoretical slope of -1/3 describing the dependence of the brush height with the salt concentration in the salted brush regime. By considering the natural counterions in the brush, the ion concentration csibrush is larger than cs, and the regime of high salt concentrations can be reached more easily. The Donnan equilibrium45 describes two solutions that are separated by a membrane. The selectively permeable membrane is constructed so that it is permeable to certain small charged ions of the solutions and impermeable to macroions. To follow this approach, we split the simulation cell into two regions (the brush and bulk regions). The boundary between these two zones is defined along the z direction by the brush thickness. We define csibulk as the free ion concentration in the bulk region. We also define csbulk ) Ns/LxLy(D - h) as the concentration of salt in the bulk region. In the case of a polyelectrolyte system, the Donnan approach is
csibulk csibrush
[ ( )]
) 1+
fNbFa
2 -1/2
(14)
2hcsbulk
The derivation of this expression can be found elsewhere.11 Equation 14 is only valid for dilute solutions. Provided that the free volume available for the free ions within the brush is reduced due to the self-volume of the polyelectrolyte chains, Kumar and Larson11 have proposed a modified version of the Donnan result
csibulk csibrush
Figure 6. (a) Brush height (h/h0) of a full charged polyelectrolyte brush as a function of the salt concentration. The solid line with a slope of (-1/3) represents the theoretical prediction for the salted brush. (b) brush Brush height dependence with respect to the (csibrush/csi,0 ) ratio, where brush csi,0 is the free ion concentration within the unsalted brush. The solid line with a slope of -0.31 represents the calculated power from the DPD results. (c) Dependence of the ratio between the free ion concentration in the bulk region and the free ion concentration within the brush as a function of the salt concentration calculated in the bulk region for three grafting densities. The dashed and solid lines represent the ratios of concentrations expected from eqs 14 and eq 15, respectively.
dependence. We may conclude that the range of salt concentrations investigated here does not correspond to that of the high salt concentrations of the brush salted regime where cs is large with respect to the counterion concentration inside of the brush. In Figure 6b, we show the dependence of the (h/h0) ratio on brush brush ), where csibrush and csi,0 (csibrush/csi,0 are the concentrations of the
)
[
( )]
fNbFa 1 1 1+ 1-m (1 - m)2 2hcsibulk
2 -1/2
(15)
At high salt concentrations, the csibulk/csibrush ratio tends to 1 in eq 14 and to 1/(1 - m) in the modified version of the Donnan equilibrium. m is an estimate of the degree of packing in the brush by considering the excluded volume of the polymer chains. In Figure 6c, we show the csibulk/csibrush ratio as a function of cbulk si calculated from DPD simulations with the curves fitted from eqs 14 and 15. Then, m is considered as an fitting parameter. The ion concentrations inside of the polyelectrolyte brush and in the bulk are given in Table 3 as a function of cs for three surface coverages. Figure 6c shows that the DPD results agree well with the modified version of the Donnan approach. At Fa ) 0.17, we see that the ratio of concentrations is larger at high salt concentrations than that expected from the standard Donnan approach. This supports the use of a parameter accounting for the excluded volume of the polymer chains. For relatively high salt concentrations, the presence of polymer beads modifies the distribution of free ions between the two regions compared to that expected from the standard Donnan approach. An increase of the ion concentration in the bulk region is observed, and this effect is much more pronounced at high grafting densities. The DPD simulations show that the concentrations calculated in the brush and bulk regions accurately reproduce the ratio predicted by the modified version of the Donnan approach. Because the high salt concentration regime is reached from salt concentrations that are greater than the polymer concentration within the
Mesoscale Modeling of Polyelectrolyte Brushes with Salt
J. Phys. Chem. B, Vol. 114, No. 21, 2010 7283 The degree of counterion condensation calculated from the integration of these curves is then reported in Table 3 and represented in Figure 7b as a function of the salt concentration. At the highest grafting density, we observe in part (a) of Figure 7 that the polyion-counterion distribution is not very sensitive to the salt concentration. We only observe a slight increase of the magnitude of the peak, whereas the shape of the distribution remains unchanged. The same distributions calculated for the smallest grafting density exhibit a significant dependence on the salt concentration. The effects of the salt concentration are much more noticeable at low grafting densities. Nevertheless, Table 3 shows that it is possible to reach degrees of counterion condensation for low-grafted brushes comparable to those obtained for high-grafted brushes; thus, the ion concentration in the bulk must slightly exceed that of the free ions in the brush. Different typical configurational of polyelectrolyte brushes are represented in Figure 8 at different salt concentrations. Conclusions
Figure 7. (a) Polyion-counterion distribution function calculated for a full charged polyelectrolyte at different grafting densities and salt concentrations. (b) Degree of condensed counterions 〈fci〉 as a function of the salt concentration at three surface coverages.
brush, it is much easier from a computational viewpoint to observe the asymptotic limit of eq 15 at small grafting densities. Figure 7a presents the polyion-counterion pair distribution functions for different grafting densities and salt concentrations.
We have performed mesoscale simulations of a polyelectrolyte brush in athermal solvent conditions in order to study the impact of the grafting density, the charge fraction, and the salt concentration on the structural properties of the brush. We have used the particle-particle particle-mesh (PPPM) method with a charge distribution function adapted to the use of a soft potential in DPD to accurately model the electrostatics. The small values of the Gouy-Chapman length explain why the counterions are confined inside of the brush for all of the surface coverages investigated here. We have also checked that the local electroneutrality is conserved along the direction normal to the grafted surface. The DPD simulations predict a parabolic shape for the monomer density profiles in contrast to that predicted by previous computer simulations where the solvent is not treated as individual particles. The parabolic profile of polyelectrolyte brushes calculated from DPD simulations is in line with that determined from neutron reflectivity measurements and density functional theory calculations. Our simulations establish an increase of the brush thickness with the
Figure 8. Typical configurations of a polyelectrolyte brush with added salt. The surface coverage Fa is fixed to 0.51, and the salt concentration cs is equal to (a) 0.014 and (b) 0.374.
7284
J. Phys. Chem. B, Vol. 114, No. 21, 2010
Ibergay et al.
grafting density. The dependence of the brush height with respect to the grating density is found to be identical to that predicted by the scaling law of the nonlinear osmotic brush regime. In the case of a neutral polymer brush, the theoretical dependence of the brush height with the surface coverage is also reproduced. We have also analyzed the impact of increasing the charge fraction from 0 to 1 in polyelectrolyte brushes. We have shown that the DPD simulations accurately reproduce the theoretical dependence of the brush height with respect to the charge fraction in the nonlinear osmotic brush regime. We have found that most of the counterions (between 97 and 100%) are confined to the brush layer, but only a specific fraction of them (measured by the degree of counterion condensation) are in close association with the polymer beads. This degree of counterion condensation is close to 99% for full charged polymer chains at high grafting densities and decreases as the grafting density is reduced. We have investigated the effect of the salt concentration of polyelectrolyte brushes formed by full charged polymer chains at three different surface coverages. Modeling the polyelectrolyte brush at a high surface coverage has allowed us to reproduce the transition from the nonlinear osmotic brush regime where h is nearly independent of cs to the salted regime where the brush height shrinks upon the addition of salt. We have also shown that it is possible to reproduce the theoretical slope of -1/3 on the condition that the salt concentration is calculated in the brush region (csibrush) and not in the whole simulation cell (cs). This indicates that we reach the regime of high salt concentrations with relatively small values of (cs). We have also calculated the distribution of the free ions in the brush and bulk regions and compared these distributions with a modified version of the Donnan equilibrium. We have demonstrated that the use of the PPPM method for the calculation of the electrostatic forces within the conventional soft DPD forces allows one to reproduce the different dependences of the brush height with respect to the grafting density, charge fraction, and salt concentration. As a result, we are currently applying the DPD methodology for the study of the interaction between two polyelectrolyte brushes as a function of the separation distance. Acknowledgment. The authors want to acknowledge R. D. Groot and F. Goujon for helpful discussions. This work was granted access to the HPC resources of IDRIS under the allocation 2009-i2009092119 made by GENCI (Grand Equipement National de Calcul Intensif). Appendix A: Electrostatic Potential (PPPM Method) The empirical expressions of the electrostatic potential used in the PPPM method21 are given by the following equations 4πreU(r) ) Γrc
{
( ) ( ) ( )
( )
rrc 2 rrc 4 52 4 rrc 2 + - 0.13587 35 5 re 5 re re rrc 6 re - 3.2100 1 rrc 2re rrc re
5.145
(r < re /rc) (re /rc < r < 2re /rc) (r < re /rc)
(A-1)
The potential and the corresponding force are represented in Figure A-1. We add for comparison the Coulombic potential and force, which diverge at r ) 0.
Figure A-1. Electrostatic potential (left axis) and force (right axis) calculated from the PPPM method. The x-axis is expressed in standard DPD units (r/rc), whereas the top axis gives the distance r/re. The potential and force expressions are plotted for two equal-sign charge distributions.
The depth of the potential for two opposite-sign charge distributions is equal to (52/35)(Γrc)/(4πre) ) 1.64rc/re. Two ions form an ion pair if the interaction energy between them is much lower than kBT. This implies that the ion pairing might occur if 1.64rc/re < 1 or re < 1.64rc. The two ions dissociates if re > 1.64rc. It means, that the value of re ) 1.64rc allows a short-range ion-ion repulsion in the DPD model that can counteract the ion pairing. Appendix B: Thermal and Mechanical Equilibria The method of planes (MOP)36 introduced by Todd, Evans, and Daivis is designed to calculate average cross-sectional pressures. The total pressure is the sum of the kinetic pkin and potential pconf contributions
pRz(z) )
〈
〉
N N mivR,i sgn(Vi,z) 1 1 + f sgn(zi - z)〉 〈 A i)1 δt 2A i)1 R,i (B-1)
∑
∑
where vR,i is the R component of the velocity of particle i and fR,i is the R component of the total force on particle i. The first term of eq B-1 represents the kinetic part which is due to the momentum of the molecules as they cross the area during ∆t. If between time t and t + ∆t particle i moves through planes, we use the sign of the z component of the velocity to specify the direction of the crossing. The second term represents the configurational part of the pressure and sums the conservative and electrostatic contributions. Figure B-1 shows the normal component of the kinetic term along the z direction in a polymer brush. We observe that the profile of the kinetic part is constant at each z, as expected for a system at thermal equilibrium. This profile leads to an average temperature of 1.05, which agrees well with the input temperature of 1.0. Panel (b) of this figure shows that the normal component of the pressure shows a flat profile along the direction normal to the surface. We also check that the electrostatic contributions (right axis) compensate exactly for
Mesoscale Modeling of Polyelectrolyte Brushes with Salt
J. Phys. Chem. B, Vol. 114, No. 21, 2010 7285
Figure B-1. Normal component of the (a) kinetic and (b) configurational parts of the pressure calculated from the method of planes in a polymer brush with cs ) 0.239 and Fa ) 1.01. The total part pkin + pC + pele given in (b) (left axis) is decomposed into the sum of conservative pC, kinetic pkin, and electrostatic pele contributions. The pC contribution is calculated with the conservative force fijC defined in eq 2.
the profile of pkin(z) + pconf(z) to ensure the mechanical equilibrium of the polymer brush. References and Notes (1) Misra, S.; Varanasi, S.; Varanasi, P. P. Macromolecules 1989, 22, 4173. (2) Pincus, P. Macromolecules 1991, 24, 2912. (3) Borisov, O. V.; Zhulina, E. B.; Birshtein, T. M. Macromolecules 1994, 27, 4795. (4) Misra, S.; Tirrell, M. Macromolecules 1996, 29, 6056. (5) Naji, A.; Netz, R. R.; Seidel, C. Eur. Phys. J. E 2003, 12, 223. (6) Ahrens, H.; Fo¨rster, S.; Helm, C. A.; Kumar, N. A.; Naji, A.; Netz, R. R.; Seidel, C. J. Phys. Chem. B 2004, 108, 16870. (7) Chen, H.; Zajac, R.; Chakrabarti, A. J. Chem. Phys. 1996, 104, 1579. (8) Csajka, F. S.; Seidel, C. Macromolecules 2000, 33, 2728. (9) Seidel, C. Macromolecules 2003, 36, 2536. (10) Hehmeyer, O. J.; Stevens, M. J. J. Chem. Phys. 2005, 122, 134909. (11) Kumar, N. A.; Seidel, C. Macromolecules 2005, 38, 9341. (12) Hehmeyer, O. J.; Arya, G.; Panagiotopoulos, A. Z. J. Chem. Phys. 2007, 126, 244902. (13) Balastre, M.; Li, F.; Schorr, P.; Yang, J.; Mays, J. W.; Tirrell, M. V. Macromolecules 2002, 35, 9480. (14) Romet-Lemonne, G.; Daillant, J.; Geunoun, P.; Yang, J.; Mays, J. W. Phys. ReV. Lett. 2004, 93, 148301. (15) Argillier, J. F.; Tirrell, M. Theor. Chim. Acta 1992, 82, 343. (16) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Europhys. Lett. 1992, 19, 155. (17) Koelman, J. M. V. A.; Hoogerbrugge, P. J. Europhys. Lett. 1993, 21, 363. (18) Goujon, F.; Malfreyt, P. J.; Tildesley, D. J. Chem. Phys. 2008, 129, 034902. (19) Lahmar, F.; Tzoumanekas, C.; Theodorou, D. N.; Rousseau, B. Macromolecules 2009, 42, 7485.
(20) Goujon, F.; Malfreyt, P.; Tildesley, D. J. Macromolecules 2009, 42, 4340. (21) Groot, R. D. J. Chem. Phys. 2003, 118, 11265. (22) Gonzalez-Melchor, M.; Mayoral, E.; Velazquez, M. E.; Alejandre, J. J. Chem. Phys. 2006, 125, 224107/1. (23) Ibergay, C.; Malfreyt, P.; Tildesley, D. J. J. Chem. Theory. Comput. 2009, 5, 3245. (24) Malfreyt, P.; Tildesley, D. J. Langmuir 2000, 16, 4732. (25) Goujon, F.; Malfreyt, P.; Tildesley, D. J. ChemPhysChem 2004, 5, 100. (26) Irfachsyad, D.; Tildesley, D. J.; Malfreyt, P. Phys. Chem. Chem. Phys. 2002, 4, 3008. (27) Goujon, F.; Malfreyt, P.; Tildesley, D. J. Mol. Phys. 2005, 103, 2675. (28) Ewald, P. P. Ann. Phys. 1921, 64, 253. (29) Yan, L. T.; Zhang, X. Soft Matter 2009, 5, 2101. (30) Yan, L. T.; Zhang, X. Langmuir 2009, 25, 3808. (31) Sirchabesan, M.; Giasson, S. Langmuir 2007, 23, 9713. (32) Espanol, P.; Warren, P. B. Europhys. Lett. 1995, 30, 191. (33) Beckers, J. V. L.; Lowe, C. P.; de Leeuw, S. Mol. Simul. 1998, 20, 368. (34) Sagui, C.; Darden, T. J. Chem. Phys. 2001, 114, 6578. (35) Groot, R. D.; Rabone, K. L. Biophys. J. 2001, 81, 725. (36) Todd, B. D.; Evans, D. J.; Daivis, P. J. Phys. ReV. E 1995, 52, 1627. (37) Irving, J. H.; Kirkwood, J. G. J. Chem. Phys. 1950, 18, 817. (38) Alexander, S. J. Phys. (Paris) 1977, 38, 983. (39) de Gennes, P. G. Macromolecules 1980, 13, 1069. (40) Jiang, T.; Wu, J. J. Phys. Chem. B 2008, 112, 7713. (41) Naji, A.; Seidel, C.; Netz, R. R. AdV. Polym. Sci. 2006, 198, 149. (42) Milner, S. T. Science 1991, 251, 905. (43) Pal, S.; Seidel, C. Macromol. Theory Simul. 2006, 15, 668. (44) Zhulina, E. B.; Wolterink, J. K.; Borisov, O. B. Macromolecules 2000, 33, 4945. (45) Donnan, F. G. Z. Elektrochem. 1911, 17, 572.
JP9115832
