9802
J. Phys. Chem. B 2008, 112, 9802–9809
Simulating Equilibrium Surface Forces in Polymer Solutions Using a Canonical Grid Method Martin Turesson,*,† Clifford E. Woodward,‡ Torbjo¨rn Åkesson,† and Jan Forsman† Theoretical Chemistry, Chemical Center, P.O. Box 124, S-221 00 Lund, Sweden, and UniVersity College, ADFA, Canberra ACT 2600, Australia ReceiVed: March 8, 2008; ReVised Manuscript ReceiVed: May 8, 2008
A new simulation method for nonuniform polymer solutions between planar surfaces at full chemical equilibrium is described. The technique uses a grid of points in a two-dimensional thermodynamic space, labeled by surface area and surface separations. Free energy differences between these points are determined via Bennett’s optimized rates method in the canonical ensemble. Subsequently, loci of constant chemical potential are determined within the grid via simple numerical interpolation. In this way, a series of free energy versus separation curves are determined for a number of different chemical potentials. The method is applied to the case of hard sphere polymers between attractive surfaces, and its veracity is confirmed via comparisons with established alternative simulation techniques, namely, the grand canonical ensemble and isotension ensemble methods. The former method is shown to fail when the degree of polymerization is too large. An interesting interplay between repulsive steric interactions and attractive bridging forces occurs as the surface attraction and bulk monomer density are varied. This behavior is further explored using polymer density functional theory, which is shown to be in good agreement with the simulations. Our results are also discussed in light of recent self-consistent field calculations which correct the original deGennes results for infinitely long polymers. In particular, we look at the role of chain ends by investigating the behavior of ring polymers. 1. Introduction An interesting property of polymers is their ability to regulate surface forces and colloid stability.1,2 The net effect of polymer addition depends on a number of parameters, such as average molecular weight, polydispersity, chain rigidity, solvent quality, and the effective interactions between the monomers and the surface of the colloid particles. There were several early attempts to theoretically rationalize polymer behavior,3–8 but it has been subsequently shown that, in many cases, the theoretical models, or approximations contained in them, were too simplistic. The treatment by deGennes,8,9 utilizing the Edwards theory, remains popular, despite its approximate nature. In particular, it relies upon the infinite chain limit. This assumption may lead to predictions that are even qualitatively incorrect,10 including a theorem stating that8 polymer-induced surface interactions are always attractive at full equilibrium. The Edwards-deGennes theory has been extended11–13 to include the effects of finite chain length. One of the predictions of this corrected theory is that dissolved polymers may generate a repulsive surface force even at full equilibrium and that this repulsion is caused by the presence of chain ends. Given that these theories are approximate in nature, it is desirable to investigate polymer-mediated surface forces by more accurate methods. Simulations are, in principle, the best choice for investigating a given model system, as (statistical noise notwithstanding) they provide an exact solution. Unfortunately, simulations are also computationally very demanding, which limits the types of model systems that can be investigated in practice. In this work, we will present a new method for * Corresponding author. † Chemical Center. ‡ University College.
simulating equilibrium surface forces, which is particularly useful for polymer models. It is a general method, that does not require molecule insertions or deletions in order to establish bulk equilibrium. Our method is closely related to techniques utilizing the isotension ensemble14–17 but offers some advantages. For example, it involves many independent simulations and is thus particularly advantageous if computer clusters are used. However, the scope of this work is not just limited to the introduction of a new simulation method. We apply the method to calculate the interaction between adsorbing surfaces in polymer solutions. Our results demonstrate an interesting interplay between surface-monomer adsorption strength and polymer concentration, as they affect these surface interactions. We will also confirm that density functional theory due to Woodward,18 combined with a generalized Flory-dimer treatment of excluded volume interactions, generates remarkably accurate predictions for surface forces in the systems under investigation. In the final part of our work, we will show that the new method is also well-suited for models where the chain architecture is more complex than the simple linear form. In particular, we will investigate ring polymers, which are particularly interesting as they contain no end monomers. Thus, we are able to test the assertion that repulsive corrections to the deGennes approach is entirely due to chain ends. 2. Model We model two interacting colloidal particles19 via the Derjaguin approximation;20 that is, they are treated as two infinite, flat, and parallel surfaces. Periodic boundary conditions apply in the lateral (x, y) directions, and the simulation box has
10.1021/jp8020529 CCC: $40.75 2008 American Chemical Society Published on Web 07/18/2008
Equilibrium Surface Forces in Polymer Solutions
J. Phys. Chem. B, Vol. 112, No. 32, 2008 9803
Figure 1. (a) The polymers are modeled as pearl-necklace chains, with a diameter, σ. Solutions containing either linear or ring polymers have been simulated. The chain configurations of a fixed number of polymers are allowed to relax between two planar walls separated with a distance h. The exchange with the bulk is not explicitly performed but the polymer chemical potential is nevertheless kept constant as h changes; i.e., µ(slit) - µ(bulk) ) 0. (b) Two types of walls have been modeled. One we refer to as strongly adsorbing (Aw ) 50), and the other we refer to as moderately adsorbing (Aw ) 25). The strength is regulated by the value of Aw.
a side length, L, and a surface area, S ) L2. As long as the surface separation h and the “size” of the polymers are small compared with the diameter of the colloidal particles, this is expected to be a reasonable model system. The polymers are modeled as pearl-necklace chains, constituted of hard sphere monomers with a diameter, σ, which also defines the constant bond length (see Figure 1). A McMillan-Mayer approach is adopted, where the (good) solvent is implicit. The effective pair potential, U(R), between two monomers, separated by a distance R is thus described by:
U(R) )
{
0 , R>σ ∞ , Reσ
(1)
The internal chain configurations are efficiently sampled via reptation and crank-shaft rotations. In the latter case, a monomer is randomly chosen, and an attempt is made to rotate it around an axis formed by its adjacent monomers on each side. In an accepted reptation step, a randomly chosen end monomer is moved to the other end of the chain. We also included global translations of whole polymers. In this study, we simulate interactions between particles with polymer-adsorbing surfaces. Specifically, all monomers are attracted to the surfaces via a truncated and shifted adsorption potential:
V(z, h) ) w(z) + w(h - z) where,
βw(z) )
[
[( σz ) - ( σz ) ] + w
Awe-z⁄σ 0
6
3
0
(2)
z e zc z > zc
(3)
We define z as the distance from the left surface and zc ) 2σ. The constant w0 is chosen such that w(zc) ) 0. The strength of the adsorption potential is given by the value of Aw. In this study, we have considered two choices: Aw ) 50, providing “strong” adsorption, and Aw ) 25, which we denote as a “moderate” adsorption potential. These adsorption potentials are depicted in Figure 1b. We aim to model an open system, allowing diffusive equilibrium between the bulk and slit regions to be maintained, as the separation is changed. Our goal is to compute the equilibrium density, corresponding to a constant polymer
Figure 2. Each point, labeled [Sl; hm], in the depicted grid, corresponds to a canonical simulation with a density determined by the values of the surface area, Sl, and the separation, hm. Starting from a reference point, the locus of constant polymer chemical potential (curved line) is determined.
chemical potential, µ, as a function of the surface separation. This kind of boundary condition has been treated using Monte Carlo methods in the grand canonical ensemble,21 and the isotension ensemble (IE). In the latter case, free energy changes have been calculated using perturbation theory15 and an optimized rates method.16 These methods have already been implemented for solutions containing neutral22,23 as well as charged polymers.24 Here, we will consider an alternative approach, where a number of canonical ensemble simulations are performed for varying surface area and separation. Within this large grid of data points, free energy differences, and the equilibrium pressure acting parallel with the surfaces, P|, can be extracted to construct loci of constant chemical potential. Thus, there is no need for polymer insertions or deletions, which makes the method more efficient than grand canonical techniques (even when configurational bias techniques25 are employed). We shall call this approach the canonical grid method (CGM); it is described in the next section. 3. Canonical Grid Method In the CGM, we generate a number (NS × Nh) of canonical simulations, each containing the same number of polymers, Np, with a degree of polymerization denoted r. NS is the number of simulated surface areas. For each area, Nh surface-surface separations are simulated. The systems span a rectangular grid where each point corresponds to a certain density, as specified by the surface area and surface-surface separation; see Figure 2. Subsequent to the simulations, free energy differences between all neighboring grid points are evaluated using Bennett’s free energy variance minimization method.26 By interpolating (using a cubic spline function) between the points, free energy curves corresponding to different bulk densities are obtained. The first issue to deal with is the size of the grid. We start by choosing a set of areas: nS ) {S0, S0fS,..., S0fSNs-1}. For each of these areas, we simulate a set of separations distributed as nh ) {h0, h0fh,..., h0fhNh-1}. The factors fSk-1 and fhm-1 give the relative change in size of the grid point labeled (km), compared with the smallest simulated area (S0) and separation (h0). 3.1. Bennett’s Method. This is an optimized method to determine free energy differences between two overlapping canonical ensembles.26,27 By calculating the ratio between the
9804 J. Phys. Chem. B, Vol. 112, No. 32, 2008
Turesson et al.
Figure 3. Schematic picture showing perturbations (thick arrows) in S and h, evaluated at an arbitrary point in the grid, labeled (km). During the simulation at this point, 〈∆U(Sk f Sk(1)〉 and 〈∆U(hm f hm(1)〉 are determined. Subsequent to all grid point simulations, corresponding averages are utilized in Bennett’s method to compute accurate free energy differences within the grid.
two partition functions with help of a weight function, Bennett’s method minimizes the variance of the free energy difference between the systems. In all grid points, ensemble averaged perturbation energies, 〈∆U(Si f Si(1)〉ij and 〈∆U(hj f hj(1)〉ij, are evaluated by virtual scaling of the particle coordinates parallel and transverse to the alignment of the surfaces. For a detailed description, see ref 23. As an example, consider a perturbation (hm f hm+1) in the point labeled (km); see Figure 3. Bennett’s expression for the quotient between the partition functions Qm and Qm+1 is given by
[
]
Qm+1 〈F(∆Um + C)〉m,k ) eC Qm 〈F(∆Um+1 - C)〉m+1,k
(4)
where ∆Um ) Um+1 - Um and ∆Um+1 ) Um - Um+1. Furthermore, F is the Fermi function, F(x) ) 1/(1 + ex). By virtual transitions from hm to hm+1, we map the system at hm+1 and vice versa. Equation 4 is valid for any value of the constant C, but the value that minimizes the variance in the Helmholtz free energy difference is given by
C ) ln[Qm+1/Qm] ) -β∆Amfm+1
Figure 4. Constant µ locus in Figure 2 is found via two consecutive thermodynamic paths, denoted I and II. (a) Process I is performed at constant area (Sref). The free energy difference is ∆A(I) ) aref S (hj) aref S (href). (b) Process II is performed at constant separation (hj) resulting in ∆A(II) ) ajh(Seq) - ajh(Sref). A cubic spline function is used to interpolate between the data points.
convenient to choose one of the nS areas that have been simulated. The choice of reference area determines the bulk density Fm. Without loss of generality, we shall define auxiliary free energy differences, akS(h) ≡ A(Sk, h) - A(Sk, h0) and akh(S) ≡ A(S, hk) - A(S0, hk). These are depicted in Figure 4. As we change the separation from href to some arbitrary separation hj, we wish to keep µ constant:
Np∆µ ) ∆A + ∆(P|Sh) ) 0
(6)
Here, P| is the average value of the pressure tensor component acting parallel with the surfaces.29 Furthermore, ∆A ) A(hj, S) - A(href, Sref), with a similar expression for the change in the pressure term. This criterion is fulfilled for one specific “equilibrium” surface, denoted Seq, at hj. Let us calculate ∆A as a two-step process (I + II); see Figure 4. The first step (I, Figure 4a) involves the free energy difference between href and hj, at constant area (Sref), whereas the second step (II, Figure 4b) is the difference in free energy, between Sref and Seq at constant separation (hj). From the simulations, ah(S) for all surface separations and aS(h) for all surface areas are obtained, using eq 4. Note that ah(S) and aS(h) are constructed by adding free energy differences obtained from perturbations between adjacent points in the grid. The relevant free energy differences are
(5)
Note that Bennett’s method requires a two-point evaluation; that is, averages at both initial and final states are needed in order to obtain ∆A. Thus, subsequently to all grid simulations, eq 4 is solved iteratively in a self-consistent way, such that C ) ln[Qm+1/Qm]. By using this technique repeatedly (for both h and S perturbations), optimized free energy differences between all points in the grid are obtained. Bennett’s method has also been coupled with isotension ensemble simulations, in order to establish equilibrium surface forces.16,22–24,28 In order to span a reasonable range of surface separations with a simple rectangular grid, we have to run quite a few canonical simulations. On the other hand, each of these “grid points” will usually only require a rather short simulation time. In other words, our approach is inherently well-adapted to parallel runs on a cluster (“trivial parallelisation”). 3.2. Finding the Constant µ Locus. When the points of the grid have been simulated, we are in a position to construct surface force curves. We start by choosing a reference separation, href from which we construct all free energy differences. Then, we choose a reference surface area, Sref. In principle, this can be any area in the range between S0 and S0fNS S-1, but it is
This division of ∆A into subprocesses is of course not unique, and there are many other options (e.g., the reverse process, II + I). This suggests the possibility of statistical improvement to the estimated constant µ locus, an aspect we have not explored in this introductory study. However, we have checked that results obtained for different choices of reference separations agree. The parallel pressure, P|(hi, S), at a separation hi and a surface area S is given by
( )
i 1 ∂ah(S) P|(hi, S) ) hi ∂S
Np
(8)
It is conveniently obtained by numerical derivation of ah(S) with respect to S. In order to accurately determine Seq at hj, by use of eq 6, interpolation between simulated data points in ajh(S) and P|(hj, S) is made by a cubic spline function.30 By successively increasing the separation, relative to the reference separation, eq 6 is solved for Seq, and the locus of constant polymer chemical potential is determined.
Equilibrium Surface Forces in Polymer Solutions
J. Phys. Chem. B, Vol. 112, No. 32, 2008 9805
Figure 5. (a) Set of osmotic pressure curves, ranging from Fmσ3 ≈ 0.017 to Fmσ3 ≈ 0.13, obtained with the CGM. Two curves obtained with the IE method, denoted by symbols, are included for comparison. The walls are moderately adsorbing, Aw ) 25, and the solutions contain 20-mers (r ) 20). (b) 3D surface of the corresponding net osmotic pressures, βPnetσ3.
4. Calculating Surface Interactions We report surface interactions either in terms of the osmotic pressure, P, or the free energy per unit area, gs. The bulk density is given by the choice of reference point, Sref, href, within the grid. The osmotic pressure, acting perpendicularly to the surfaces, can be obtained as a grand potential derivative:
P)
d(P|(h)h) dh
(9)
An important consistency (equilibrium) check is that the same quantity is obtained via the “virial” expression:
P)
∫0h Fm(z)( dw(hdh- z) ) dz
(10)
where Fm(z) is the average monomer density at a position z in the slit. This is a sensitive test, and even minor errors will generally lead to a measurable discrepancy between eq 9 and eq 10. The interaction free energy per unit area is calculated as
gs ) (Pb - P|(h))h
(11)
where Pb is the osmotic pressure in the bulk, estimated from P at large surface separations. Moreover, in this regime, gs will approach twice the surface tension at a single surface, γs. Hence, we will report the net free energy per unit area as ∆gs(h) ≡ gs(h) - 2γs. 5. Comparison with Other Techniques 5.1. Isotension Ensemble. The CGM is closely related to a transition rates method utilizing the isotension ensemble, IE.16,22,23 The IE approach relies upon probability distributions expressed in terms of the area fluctuations generated in the simulations. These fluctuations are balanced by an externally applied parallel pressure. The equilibrium parallel pressure at a different surface separation is obtained via virtual perturbations of the separation. Free energy differences are conveniently computed with an IE adapted version of Bennett’s method.16 This method has been successfully used to produce equilibrium surface interactions in solutions containing simple molecules,16,28 as well as neutral22,23 or charged polymers.24 CGM offers some potential advantages, in comparison with this technique. IE simulations can be problematic in the vicinity of phase boundaries. Phase transitions can be suppressed by switching to a “restricted” IE,31 albeit at the cost of greater computational complexity and some loss of statistical efficiency. In the present work, we study hard sphere systems, with no possibility of fluid phase transitions. The accuracy was found to be about as good
as that with the IE technique. On the other hand, from one grid simulation, a whole set of osmotic pressure curves is obtained, ranging at least 1 order of magnitude in terms of bulk density. 5.2. Grand Canonical Ensemble. Grand canonical ensemble simulations (GCE) are well-suited for studies of open systems but are, in practice, limited to rather simple molecules. When particles are linked, forming long polymers, insertions and deletions are thwarted by very low acceptance probabilities. This can, to some extent, be circumvented by the use of biased insertion/deletion techniques.23,32–34 Still, even these methods fail to work at fairly modest polymer lengths, as we will demonstrate. Furthermore, complex chain architectures, such as ring polymers, are problematic to handle in a GCE simulation (especially with rigid bonds). By contrast, ring polymers are relatively easily incorporated into the IE and CGM approaches. 6. Results In this section, we report on surface interactions computed with the CGM and the IE technique. Some comparisons will include GCE results. The surfaces attract monomers, as described by eq 3. We will consider surface interactions mediated by polymers with r ) 20 and r ) 50, in a slit between moderately (Aw ) 25) and strongly (Aw ) 50) adsorbing surfaces. The number of simulated grid points for each system exceeds 2000 (NS ≈ 35 and Nh ≈ 70). All grids were constructed using a “reference separation”, href ) 6.04 σ. The specific choice of Sref then defines the bulk density and the corresponding surface force curve. Furthermore, Np ) 16. In most cases, data points are more frequent than symbols in the graphs. 6.1. Moderately Adsorbing Surfaces. The distributions of areas and separations were generated, using fS ) 1.03 and fh ) 1.016. Furthermore, we chose S0 ) 275 σ2 and h0 ) 4 σ. By choosing different reference areas in the grid, the osmotic pressure curves in Figure 5 were generated, via eq 9. The symbols denote results obtained with the IE approach. The bulk monomer density, Fmσ3, ranges from 0.017 to 0.13. For the highest concentration, the size of the grid is too small to cover separations larger than 7 σ (it is, however, possible to subsequently extend the grid to increase the range). The agreement between results obtained from the CGM and the wellestablished IE approach confirms the usefulness and accuracy of the new method. In Figure 5b, a three-dimensional (3D) plot of the corresponding net osmotic pressures, Pnet ) P - Pb, is shown. All systems display bridging-dominated attractions.35 6.2. Strongly Adsorbing Surfaces. If the monomers adsorb strongly to the surfaces, the picture changes quite dramatically. For low bulk monomer densities, we still observe net attractive osmotic pressures. For higher densities, however, the interactions
9806 J. Phys. Chem. B, Vol. 112, No. 32, 2008
Turesson et al.
Figure 6. (a) Set of osmotic pressure curves, for strongly adsorbing surfaces (Aw ) 50), as obtained with the CGM. The solution contains 20-mers, and the bulk monomer concentration ranges from Fmσ3 ≈ 0.003 to Fmσ3 ≈ 0.15. Results from IE simulations, for two different bulk densities, are given by symbols. (b) 3D surface plot (in the range from 5-6 σ) of the corresponding net osmotic pressures. (c) Net interaction free energies ∆gs, for various bulk concentrations, in the transition regime where the interaction turns repulsive. Isotension ensemble results at intermediate and high bulk densities are also included for comparison.
turn strongly repulsive. The effect is illustrated in Figure 6, showing results from a grid simulation of 20-mers and Aw ) 50. In this case, we chose fh ) 1.016, fS ) 1.012, h0 ) 4 σ, and S0 ) 200 σ2. The bulk monomer density ranges from Fmσ3 ≈ 0.003 to Fmσ3 ≈ 0.15. Our interpretation of the observed response is that, at high concentrations and with strongly adsorbing surfaces, the surfaces become saturated by polymers. Bridging chains will then find few “available sites” at the opposing surface. Hence, the repulsive effects on an increased density in the slit as the surfaces approach, that is, excluded volume and ideal entropy, will dominate the attractive bridging term. These results contradict a well-known theorem by deGennes, predicting that all polymer mediated surface interactions are purely attractive under full equilibrium conditions.8 The overall behavior is illustrated in a 3D plot in graph b of Figure 6. In Figure 6c, we focus on the net interaction free energy, in the concentration regime where the interaction turns repulsive (upon increasing Fm). We have also included an example where the bulk density is very high. In this case, the repulsion is actually weaker at large separation, as compared with a bulk with intermediate concentration. The reason is that, above the overlap concentration, the range of the surface interaction is limited by the correlation length, which decreases as the concentration increases. In other words, at high concentrations, the repulsion is strong but short-ranged. We shall highlight this nonmonotonic behavior in more detail below, using density functional calculations. In Figure 7, we illustrate that optimized noninsertion methods can be used to calculate equilibrium surface interactions, even when the walls are strongly adsorbing and the chains are relatively long. The IE method is used to generate the net interaction free energy in a 50-mer solution. Longer chains adsorb more strongly, thereby producing rather dense layers even at low concentrations. This leads to an excluded volume repulsion between nearly saturated surfaces. This effect will be further investigated by DFT calculations below.
Figure 7. Net interaction free energy between two strongly adsorbing (Aw ) 50) surfaces in contact with a bulk solution containing 50-mers, at a monomer concentration of Fmσ3 ≈ 0.01. The data were obtained via IE simulations.
6.3. Breakdown of the GCE Approach. In the configurationally biased GCE approach, polymers are gradually inserted, according to an algorithm that is biased by a Boltzmann weight. This bias is corrected for at the end of the insertion/deletion step.32–34 The technique can be used to construct surface free energies in polymer solutions, at least if the chains are linear, flexible, and relatively short, and with surfaces that are not too strongly adsorbing.23,36 In Figure 8, we compare the IE method and the grand canonical simulations for a system containing 50-mers at Aw ) 50. The strong repulsion observed with the IE method is not captured with the GCE approach. Note that the agreement between the osmotic pressure, as obtained as a free energy derivative, eq 9, and that from the virial route, eq 10, ensures that the IE simulations provide the correct equilibrium surface interaction. Furthermore, no hysteresis was found. An attempt to increase the system size from 〈Np〉 ≈ 20, to 〈Np〉 ≈ 80, was made in the GCE. Still no repulsive forces could be detected. Previous studies have demonstrated that GCE simulations are even more problematic when the chains are stiff. In such cases,
Equilibrium Surface Forces in Polymer Solutions
Figure 8. Osmotic pressures calculated with IE simulations and configurationally biased GCE simulations. With moderate adsorption potential and short polymers (see inset), the agreement between the methods is excellent. When the adsorption strength and the chain length are increased to Aw ) 50 and r ) 50, the GCE approach breaks down; i.e., the repulsion at short separations is not captured at all. The IE data has been obtained by using eqs 9 and 10 and stepping both upward (from h ) 4 Å) and downward. Fmσ3 ≈ 0.01.
the “breakdown” occurs already for shorter chains (even if they are ideal), and with weaker adsorption potentials.23 It might be possible to establish correct equilibrium densities for longer chains by the use of grand canonical expanded ensemble methods, where a chain is allowed to grow or diminish monomer-by-monomer during the simulation.37,38 However, the existence of a partial chain in the system will lead to incorrect thermodynamic quantities (an error that diminishes with increasing system size). Hence, one would then either have to restrict pressure calculations to configurations that only contain complete chains, or else establish the correct densities and then
J. Phys. Chem. B, Vol. 112, No. 32, 2008 9807 calculate pressures in separate canonical simulations. The first approach will lead to poor statistics and the latter is cumbersome. 6.4. Density Functional Calculations. In this section, we will compare our simulation results with predictions by polymer density functional theory, DFT. Specifically, we will adopt the formulation presented by Woodward.18 This theory is exact for ideal chains, with point-like monomers. In other words, all configurations are accounted for, subject to a mean field. This is an important property of the theory, as we shall see below. In our model system, the monomers are treated as hard spheres, and excluded volume effects have to be approximated. Previous simulation comparisons have demonstrated that a version based on the generalized Flory-Dimer treatment23,39–42 is accurate, and we shall adopt it in this work. Since the Woodward theory is quite standard and has been thoroughly described elsewhere40,42,43 (with some minor differences of the way in which the hard cores are treated), we will not present the formal aspects here but rather focus on the results of our calculations. The DFT predictions of interactions between strongly adsorbing (Aw ) 50) surfaces, in 20-mer solutions with various bulk densities, are presented in Figure 9. The results in Figure 9a and Figure 6a are directly comparable, and the agreement is quite satisfactory. An advantage is that the results are noisefree and that the bulk conditions are known exactly. This has allowed us to construct smooth 3D plots of the net interaction free energy, and its response to changes of the bulk concentration. These are given in Figure 9b,c. The former graph illustrates that the overall response to an increased polymer density is a more repulsive surface interaction. Still, as we mentioned earlier, there is a “saturation effect” at large separations and high densities. In this regime, a further increase of the bulk density will actually reduce the repulsion, since the bulk correlation
Figure 9. Density functional predictions of interactions between strongly adsorbing surfaces, immersed in solutions containing 20-mers of various concentrations. (a) Osmotic pressure curves. These can be compared with the corresponding simulation results, presented in Figure 6a. (b) 3D plot of net surface free energies, in the separation regime h ) 4-8 σ. This illustrates the overall repulsive response to an increased bulk concentration. (c) 3D plot of net surface free energies, in the separation regime h ) 5-6 σ. This highlights the nonmonotonic response observed at long-range, where the repulsion at high concentrations drops, upon a further increase of the bulk density.
9808 J. Phys. Chem. B, Vol. 112, No. 32, 2008
Turesson et al.
Figure 10. Comparing how surface interactions respond to changes of the bulk concentration, in solutions containing 20-mers (a) and 200-mers (b). The net surface free energies are obtained by DFT calculations. The inset in each graph is an enlargement of the response in the large separation regime.
Figure 11. Connecting the free chain ends in linear polymers, graph (a), forming ring polymers, graph (b), does not dampen the repulsive barrier. On the contrary, the onset of repulsion seems to occur at lower densities for rings compared with linear chains. The range of repulsion is diminished with ring polymers, because of the smaller radius of gyration. r ) 20, Aw ) 50.
length drops, generating a strong but more short-ranged repulsion. We have discussed this in connection to Figure 6c, but the effect is more clearly illustrated in Figure 9c. Another advantage of the DFT approach is that the calculations run many orders of magnitude faster than corresponding simulations. This permits studies of long (even infinite)10 chains, and it is of considerable interest to make corresponding analyses of solutions containing relatively long polymers. In Figure 10, we make such a comparison, between 20-mer and 200-mer solutions. On a large scale, the response to bulk density changes appears to be stronger with the shorter chains. However, this is partly misleading. In the insets, which focus on the large separation regime, we observe a qualitatively opposite behavior, with a stronger response for the long-chain solutions. Still, the overall features are retained in the presence of long chains, that is, an attraction at low concentrations, turning repulsive as the bulk density increases. We note that, even at the lowest concentration reported, the 200-mer solution actually displays a strong repulsion at short-range. Cooperativity effects lead to very strong adsorption for long chains, which is not easily diluted away. The short-ranged repulsion is a consequence of an overlap between these adsorbed layers and occurs at the monomer-monomer level. Specifically, adsorbed monomers will experience excluded volume interactions from the adsorbed layer at the opposing surface when the separation is below 5 σ. This repulsion persists down to very low concentrations. Only when Fmσ3 ≈ 1 × 10-8 or below, the interaction free energy is wholly attractive in the 200-mer solution. 6.5. Ring Polymers. Some chain architectures impose problems when a configurationally biased GCE technique is adopted. With rigid bonds, there is a complicated constraint on the
possible coordinates in the successive insertions of the monomers of a ring polymer. This is not a problem with the CGM or the IE method, since no insertions or deletions are needed. Ring polymers are of some fundamental interest, since an extended version of the Edwards-deGennes mean-field theory predicts that any repulsion in a polymer solution is caused solely by chain ends.11–13 In this treatment, monomers are partitioned as belonging to “trains”, “loops”, or “ends/tails”, where the last category is responsible for repulsive interactions. Ring polymers then form an interesting test case, as no monomers belong to “tails”. Our predictions may also have some practical relevance, since such cyclic molecules can be synthesized. In Figure 11b, we present surface free energies in a solution containing ring polymers. These were obtained by CGM simulations, that is, we have automatically sampled a range of bulk concentrations. There is no unambiguous way to compare rings and linear chains. We have chosen to simply connect the ends, that is, retain the number of monomers per chain, as we compare the two cases. We see the same qualitative trends as we have observed with linear chains; see Figure 11a. An attraction at low concentrations is replaced by repulsive interactions at higher densities. However, the repulsion starts to dominate at lower concentrations when the ends are connected to form ring structures. Furthermore, the transition is more rapid; that is, the interactions are more sensitive to changes of the bulk density. As one would anticipate, ring polymer induced interactions are somewhat more short-ranged than their linear chain correspondence. The major effect, in terms of mediated surface forces, of connecting the ends appears to be that the molecules become more compact. It seems relevant to make connections to another
Equilibrium Surface Forces in Polymer Solutions study,35 where surface forces in solutions containing linear and star polymers were compared. The work included the same density functional theory that we have adopted as a part of this work. As the number of “arms” in the stars increases, the surface interactions become more short-ranged but also more repulsive, in qualitative agreement with our findings for rings. However, in the former case, the number of chain ends increase as arms are added, at constant molecular weight. In our case, the more compact structures contain no ends. Still, the response is similar. The strong repulsions we observe in the presence of ring polymers suggest that the extended Edwards-deGennes meanfield theory11–13 is too approximate. One possible source of error is in the description of bridging interactions. These make very strong attractive contributions to the net surface force and are counteracted (and very nearly canceled) by an opposing entropic (ideal) contribution.35 In order to properly balance these two very large terms, they need to be evaluated consistently, as is done with density functional theory, where all configurations are accounted for (and weighted by a self-consistent mean-field). It should be pointed out that, while “ends” do contribute substantially to the entropic, repulsive terms, they are also able to generate strongly attractive bridging interactions. The net outcome requires detailed and balanced analyses. 7. Conclusions We have developed a new canonical ensemble simulation method to study polymer solutions confined by planar walls, in equilibrium with a bulk solution. Our method requires a grid of thermodynamic points be constructed with the free energy differences between these points calculated accurately and efficiently. For this calculation, we employed Bennett’s optimized rates method. By appropriate thermodynamic manipulation and interpolation within the grid, a series of constant chemical potential loci can be constructed, which allow us, in this case, to calculate the interaction free energy and osmotic pressure versus surface separation. However, the method is suitably general that other types of systems and boundary conditions (e.g., constant pressure) could be treated in a similar fashion. The major advantage in the current system under study (polymer solutions) is that we are able to eliminate the need for particle insertions, which greatly facilitates the simulations, as comparisons with traditional (albeit improved) grand canonical ensemble simulations have shown. We have also shown that the polymer density functional theory is able to give reasonably good quantitative agreement with the simulations, for the parameters under study here. Both approaches confirm interesting interplay between surface attraction and polymer density on the surface interactions. In particular, at the onset of repulsion (observed for sufficiently high bulk concentrations), steric interactions more than compensate the bridging forces that dominate at lower densities. Our results for ring polymers seem to cast some doubt on recent extensions of the original Edwards-deGennes, which attribute repulsive interactions solely to polymer ends. Finally, the CGM is also admirably suited to study phase transitions in uniform systems, as the constant volume simula-
J. Phys. Chem. B, Vol. 112, No. 32, 2008 9809 tions are able, in principle, to traverse weakly metastable regions of thermodynamic space more easily than constant pressure simulations. This will allow us to study such phenomena as surface wetting and layering transitions and capillary condensation phenomena in polymer systems. References and Notes (1) Evans, F. A.; Wennerstro¨m, H. The colloidal domain: where Physics, Chemistry, Biology and Technology meet; VCH Publishers: New York, 1994. (2) Israelachvili, J. N. Intermolecular and Surface Forces, 2nd ed.; Academic Press: London, 1991. (3) Vincent, B. AdV. Colloid Interface Sci. 1974, 4, 193. (4) Napper, D. H. J. Colloid Interface Sci. 1977, 58, 390. (5) Joanny, J.-F.; Leibler, L.; deGennes, P. G. J. Polym. Sci.: Polym. Phys. Ed. 1979, 17, 1073. (6) Feigin, R. I.; Napper, D. H. J. Colloid Interface Sci. 1980, 74, 567. (7) Feigin, R. I.; Napper, D. H. J. Colloid Interface Sci. 1980, 75, 525. (8) deGennes, P. G. Macromolecules 1982, 15, 492. (9) deGennes, P. G. Rep. Prog. Phys. 1969, 32, 187. (10) Woodward, C. E.; Forsman, J. Phys. ReV. E 2006, 74, 010801. (11) Bonet-Avalos, J.; Joanny, J.-F.; Johner, A.; Semenov, A. N. Eurpohys. Lett. 1996, 35, 97. (12) Semenov, A. N. J. Phys. II France 1996, 6, 1759. (13) Semenov, A. N.; Joanny, J.-F.; Johner, A.; Bonet-Avalos, J. Macromolecules 1997, 30, 1479. (14) Parinello, M.; Rahman, A. Phys. ReV. Lett. 1980, 45, 1196. (15) Svensson, B.; Woodward, C. E. J. Chem. Phys. 1994, 100, 4575. (16) Forsman, J.; Woodward, C. E. Mol. Sim. 1997, 19, 85. (17) Attard, P. Mol. Phys. 1993, 78, 943. (18) Woodward, C. E. J. Chem. Phys. 1991, 94, 3183. (19) Dickinson, E.; Euston, S. R. J. Chem. Soc., Faraday Trans. 1991, 87, 2193. (20) Derjaguin, B. V. Kolloid Zeits. 1934, 69, 155. (21) Adams, D. J. Mol. Phys. 1975, 29, 307. (22) Czezowski, A.; Woodward, C. E. Comput. Phys. Commun. 2001, 142, 117. (23) Turesson, M.; Forsman, J.; Åkesson, T. Phys. ReV. E 2007, 76, 021801. (24) Turesson, M.; Woodward, C. E.; Åkesson, T.; Forsman, J. J. Phys. Chem. B 2008, 112, 5116. (25) Siepmann, J.; Frenkel, D. Mol. Phys. 1992, 75, 59. (26) Bennett, C. H. J. Comput. Phys. 1976, 22, 245. (27) Deitrick, G. L.; Scriven, L. E.; Davis, H. T. J. Chem. Phys. 1989, 90, 2370. (28) Forsman, J.; Woodward, C. E.; Jo¨nsson, B. Langmuir 1997, 13, 5459. (29) Hill, T. L. An Introduction to Statistical Thermodynamics; Dover: New York, 1986. (30) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. In Numerical Recipes in Fortran 77: The Art of Scientific Computing; Cambridge University Press: Cambridge, 2001. (31) Forsman, J.; Woodward, C. E. Mol. Phys. 1997, 90, 637. (32) Rosenbluth, M. N.; Rosenbluth., A. W. J. Chem. Phys. 1955, 23, 356. (33) Harris, J.; Rice, S. J. Chem. Phys. 1988, 88, 1298. (34) Siepmann, J. Mol. Phys. 1990, 70, 1145. (35) Woodward, C. E.; Forsman, J. Macromolecules 2004, 37, 7034. (36) Turesson, M.; Åkesson, T.; Forsman, J. Langmuir 2007, 23, 9555. (37) dePablo, J.; Laso, M.; Suter, U. J. Chem. Phys. 1992, 96, 2395. (38) Escobedo, F. A.; Suen, J. K. C. J. Chem. Phys. 1995, 103, 2703. (39) Honnell, K. G.; Hall, C. K. J. Chem. Phys. 1991, 95, 4481. (40) Woodward, C. E.; Yethiraj, A. J. Chem. Phys. 1994, 100, 3181. (41) Wichert, J. M.; Gulati, H. S.; Hall, C. K. J. Chem. Phys. 1996, 105, 7669. (42) Forsman, J.; Woodward, C. E. J. Chem. Phys. 2004, 120, 506. (43) Yu, Y.-X.; Wu, J. J. Chem. Phys. 2002, 117, 2368.
JP8020529
