|
EQUILIBRIUM PROPERTIES OF POLYMER SOLUTIONS AT SURFACES:
MONTE CARLO SIMULATIONS
By
JASON DE JOANNIS
A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
UNIVERSITY OF FLORIDA
2000
Copyright 2000
by
Jason de Joannis
This work I dedicate to my wife Sabrina. How little it is in comparison to her dedication
to me.
ACKNOWLEDGMENTS
I would like to recognize loannis Bitsanis for introducing me to this exciting field
of research and his constantly insightful advice and guidance. My other advisor, Chang-
Won Park, has always been a great source of stability and encouragement. I have
appreciated the efforts made by Raj Rajagopalan in my personal development as a
researcher, numerous technical discussions and also in some of the administrative
requirements of this degree. I am grateful to Anthony Ladd and Sergei Obukhov for
sharing their expertise with me through several discussions. Special thanks goes to Jorge
Jimenez for an infinite number of interesting discussions on random subjects and for the
fruitful collaboration of the past year.
Also I would like to acknowledge my parents Eugene and Diane, who have given
me constant encouragement and the foundation from which to think independently. My
brother Jeff and two sisters Alexa and Erica, have also been curious and supportive of my
research. Finally, I would like to recognize Mehmet Bayram Yildirim, for helping me
with some computational work, and Adriana Vasquez-Jimenez, through whom I have
finally reached a full understanding of "the physisorbed layer".
TABLE OF CONTENTS
page
A C K N O W L E D G M E N T S .............................................................................................. iv
A B ST R A C T .............................................................................................. viii
1. IN T R O D U C T IO N ................................................................................. 1
2. ELEMENTS OF POLYMER STATISTICAL MECHANICS.............. .................... 5
P olym er Solutions ...................................... ............................................. 5
Historical Perspective .................................... ............................ ......... 6
Ideal C hains ....................................................................... ......... 9
E excluded V olum e ........................................... .... .............................. 20
The Flory-Huggins Theorem A Fundamental Equation of State ..........................27
Scaling L aw s ............... ...... ............................ ......... ... ............ 29
P hase D iagram ......... ................................................... ............................ 32
Polymer Solutions at Hard Interfaces................ ............ .....................34
Mean-Field Theories of Adsorption on a Lattice............................... ...............34
The Ground State Dominance Solution to Continuum Mean-Field Theories...........37
Scaling D description of Polym er Adsorption ........................................ ................ 38
Square-Gradient Mean-Field Free Energy Functional .........................................39
B eyond the G round State ...................... .. ................................... .............. 41
3. EXPERIMENTAL CHARACTERIZATION OF POLYMER INTERFACES .........43
M easurement of Polymer Structure at Surfaces ...................................... ............... 43
Global Properties of Physically Adsorbed Layers............................... .............. 44
Sm all A ngle N eutron Scattering.................................. ..................................... 47
N eutron R eflectivity ...................................... ..................... ............... 49
Scattering and Reflectivity Measurements on Adsorbed Polymers .........................49
M easurement of Polymer Induced Forces ........................................ .............. 50
Surface F orce A pparatu s......................................................................... ...... 5 1
SFA M measurements in Good Solvent ......................... ......... ..................... 53
A tom ic Force M icroscopy .................................................................. ... 57
AFM Measurements on Elasticity of Single Chains and Forces Between
P olym er L ayers ................................................................... .............. 59
4. MONTE CARLO SIMULATION OF POLYMER CHAINS ON A LATTICE........61
Enumeration, an Exact Theoretical Treatment for Model Oligomers..........................62
M onte C arlo M ethods ............................................................................. .......... 64
Static M onte Carlo ............. .............................................................................64
D ynam ic M onte Carlo ................................................................ .............. 66
Configurational Bias Monte Carlo .............. ........................ ..............69
Sim ulation of B ulk Properties ...................... .. ............. .................... .............. 70
Sim u nation of A d sorption ...................... .. ............. .. ............................................ 72
5. REPTATING CONFIGURATIONAL BIAS MONTE CARLO FOR LONG CHAIN
MOLECULES AND DENSE SYSTEMS ..............................................82
Introduction ..................................... ....... .... .... .. .................... 82
Proposal and Detailed Balance of Configurational Bias Extensions ........................... 85
End-Switching Configurational Bias ......................................... .............85
Persistent End-Switching Configurational Bias ...................................................88
T ests of the Proposed A lgorithm s ...................................................... ................... 91
Static P properties ......................................................... ........................ 92
Dynamic Characteristics of the Algorithms..........................................................96
Su m m ary ...................................... ................................................. . 9 9
6. A POLYMER CHAIN TRAPPED BETWEEN ATHERMAL WALLS................. 115
In tro d u ctio n ........................................................................ ............... 1 1 5
Single Chains and Interfaces ........................................................... ..............116
Single C hain A dsorption .......................................................... .............. 117
Single C hain B ridging ........................ ........................................ 118
Single Chain Confinem ent ......................................................... ......... ..... 118
Sim ulation M mechanics ........................................................... .............. 122
R esults.................................................. ........... 124
C concentration Profiles .......................................................................... ............. 125
Force of Confinement and Connection with Surface Layer Concentration .......... 126
Universal Force Concentration Constant ......................................... ................ 127
Sum m ary ............................................................... .... ...... ......... 127
7. HOMOPOLYMERPHYSISORPTION.............................. 133
Introduction ............................................................... .. ..... ......... 133
Sim ulation M mechanics ............................................................ .............. 137
R esu lts............................ .. ................ .............. 140
Qualitative Picture of the Interface............................................... .................... 140
Quantitative Comparison to Finite Chain-Length Lattice Mean-Field Theory ...... 142
S u m m ary ................................................................................... 14 6
8. ADSORPTION OF POLYMER CHAINS WITH INFINITE MOLECULAR
W E IG H T ................... ............................................................... 16 3
Introduction ............................................................... ... ..... ......... 163
Sim ulation M echanics ............................................................ .. ........... ... 164
Theoretical Predictions ...... ...... .............. ..................................... 165
Analytical Solutions........................ ........................ ....... 165
Limiting Behavior of Lattice Mean-Field Calculations .................. .............. 167
R e su lts ......................................................... 16 8
0-Solv ent ......... ... ......................................................................... 168
Atherm al Solvent .................................... ....... ........... .............. 169
Further Details of the Layer Structure ......................................... .............172
C onclu sions ............................................................... ... ..... . ..... 174
9. COMPRESSION OF AN ADSORBED POLYMER LAYER OF FIXED MASS .. 195
In tro d u ctio n ................................................................................. 19 5
Sim ulation M mechanics .......................................................................... .............. 199
Sampling the Conformations of the Chains ............ ............................... 200
Calculation of Forces from Lattice Simulations ............. ..........................201
R esu lts ............ ......... ............. ............ ...... .. .............................. . 202
Compression of a Polymer Layer with an Athermal Wall ................................... 202
Interaction Between Two Polymer Layers......................................................209
Sum m ary ................................................................ .... ...... ......... 210
10. PRO SPECTS ................. ................. ...................... .. ..... ................ 221
Depletion at Non-Adsorbing Surfaces...................................... ... ..................221
Competitive Adsorption .............. ........... .............. .......... ......... .............. 222
Ring Polymers, Grafted Polymers and Other Problems...........................................223
11. CON CLU SION S .............. .............. ..................... ......... .................. 227
Extension of Configurational Bias Monte Carlo................................................... 227
Compression of an Isolated Chain ............ ............... ..................... ... ..............227
Test of Mean-Field and Scaling Theory of Adsorption....................... ............228
Compression of an Adsorbed Polymer Layer..................... ........................229
REFEREN CE S ............ ......... ............................ .......... .. ......... ........ 232
BIO GR APH ICAL SK ETCH ........................ ..................................... .............. 242
Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
EQUILIBRIUM PROPERTIES OF POLYMER SOLUTIONS AT SURFACES:
MONTE CARLO SIMULATIONS
By
Jason de Joannis
December 2000
Chairman: Raj Rajagopalan
Major Department: Chemical Engineering
Knowledge of conformational features of polymer solutions at interfaces is a
prerequisite for understanding and controlling important technological processes such as
colloid stabilization, biopolymer adsorption, corrosion inhibition, lubrication, adhesion,
membrane separations and aggregation-induced separations.
In this exposition, some examples are shown of the unique contributions that a
carefully designed simulation study can make towards an improved interpretation of
experimental and theoretical work. One section is devoted towards the improvement of
an existing Monte Carlo algorithm to enable better tests of theoretical results. Large-scale
end-swapping configurational bias moves provide the best chance of simulating high
chain lengths and concentrations.
The main focus of this paper is on polymers adsorbing from a good solvent onto a
solid surface that has a weak (reversible) attraction on monomers. We analyze the
concentration gradient that develops when the interface is in equilibrium with a bulk
solution, and the force induced by adsorbed chains when subject to confinement.
Strong evidence is presented supporting the scaling prediction for long chains that
the loop and overall concentration decay as z-4/3 when adsorption occurs from semi-dilute
or dilute athermal solutions. Mean-field theories predict a concentration decay that is too
rapid at all chain lengths, consequently they underestimate adsorbed amount and layer
thickness. Simulations provide the only quantitative description of finite-length random
walks with excluded volume subject to adsorption conditions.
Adsorbed chains under compression exhibit lower pressures than predicted by
mean-field theory as a result of its lower adsorption capacity. Even when interpreted with
respect scaled variables, quantitative discrepancies persist in the compression of saturated
/ semidilute layers. A different comparison shows that the simulation results presented
have substantial experimental implications. The compression of two adsorbed layers
against each other is quantitatively similar to the compression of a single layer, when the
surfaces bear high adsorbed amounts.
A study also is made of the concentration gradient and confinement energy of a
single chain between athermal walls. It is found that good agreement exists with
theoretical predictions based on the "magnetic analogy" for a universal correlation
between compression force and the depletion of segments near the walls.
Finally several remarks are made concerning further research on problems
including polymer depletion at interfaces and adsorption of polymers from a bidisperse
distribution of molecular weight.
CHAPTER 1
INTRODUCTION
In this paper, the equilibrium structure and interactions of linear homopolymer
(i.e., polymer with identical repeat units) molecules at various interfaces are examined
using a lattice Monte Carlo simulation method. Chemical and physical adsorption of
polymers on solid surfaces plays an important role in many technological processes.
Among the most important are the stabilization of colloidal suspensions (e.g., foods,
paints and inks, pharmaceuticals), adhesion, lubrication, and corrosion inhibition.
Because of the wide variety of polymers available (e.g., linear, branched, block-
copolymers, random copolymers, star, tree and homopolymers) and due to the various
surface/polymer interactions that can take place, the physical characteristics of a surface
in principle can be manipulated to any extent. One impressive example of this is in the
biological compatibilization of artificial implants.
The adsorption of linear homopolymers is the most widespread method of
stabilizing colloidal systems.1 The repulsive force that occurs when two polymer-coated
particles approach each other is termed "steric repulsion" since it arises largely from the
reduced entropy on confinement of the polymer conformations. The advantage of these
systems is that no specialized polymer synthesis or suspension preparation methods are
required. More costly methods such as the grafting of polymers onto particle surfaces are
used as a last resort and the grafting of polymers to some surfaces can be difficult.
An exact theoretical description of polymer / surface systems is overwhelming. In
principle, the equilibrium characteristics of the system depend on a multitude of factors:
monomer-monomer and monomer-solvent interaction energies (including ionic, van der
Waals forces, and hydrogen bonding forces), size and shape of the monomer and solvent
molecules, the polymer's valence angles as well as its secondary and tertiary (e.g., helical
structure of DNA) structure, and the distribution of molecular weight. Further
complications arise when the surface has a corrugated crystalline structure and can have
several kinds of long and short-range interactions with the polymer solution leading to
non-uniform adsorption characteristics.
In general no attempt is made at such a detailed description. In lieu of an exact
treatment, theories attempt to identify two or three of the most important parameters of a
problem and infer the scaling of system properties with these parameters.2 For example, a
function might be proposed, for either argumentative or observational reasons, to depend
only on two important variables x and y, f (x, y). Then the objective is to find the power
law with which these scale: f ~ x"m yn. The scaling is deduced from dimensional or
phenomenological arguments based on the relevant length and energy scales in the
system. This kind of reasoning can be surprisingly successful under asymptotic
conditions and has shed light on an amazing variety of problems in polymer physics.
When such a complex system reduces to simple functions of only a few primary variables
it is referred to as universal. That is, the results hold over a coarse grained scale
regardless of the chemical constituents of the polymer or solvent. For polymer solutions,
the chain length, N, the solvent quality, X, and the volume fraction, 0, capture most of the
information about thermodynamic properties. A discussion of theories that describe the
universal behavior of polymers at interfaces and in solution will be made in Chapter 2.
In Chapter 3 the experimental tools and observations that are relevant to the
structure and forces exhibited by polymer layers are reviewed. The most fruitful studies
have used surface force apparatus3 and small angle neutron scattering4 techniques.
Neutron scattering has been used to measure the concentration gradient within adsorbed
polymer layers. Unfortunately the limitations on this technique and the lack of atomistic
modeling capacity prevent simulations from being quantitatively compared with
corresponding measurements. Instead simulations are used to test the coarse-grained
theories of polymer adsorption which have remained unproven for two decades (see
Chapters 7 and 8).
Surface force apparatus measurements have succeeded in measuring the forces
between surfaces bearing adsorbed polymer chains. An understanding of the detailed
structure of the interface in terms of its concentration gradient and conformational
architecture is a prerequisite for making correct interpretations of these measurements. In
Chapter 9 we investigate the force vs. distance and the corresponding conformational
changes that take place. We make an interpretation based on the properties of layers of
unperturbed layers of Chapter 7 and 8. This is illustrative of the fact that simulations can
be used as a great interpretive tool.
In Chapter 6 we discuss a model problem that has found surprising application.
Theories and simulations of single chains abound in literature because they can be
studied and understood more rigorously. Such is the case with the compression of a
single chain between non-adsorbing walls. The lateral expansion of the chain has been
described by scaling and other theoretical arguments. This may appear to be a problem of
little consequence other than academic. However the lateral expansion of an adsorbed
chain into a pancake is subject to precisely the same physical arguments.2 The confining
surface is substituted by the strength of the adsorption field.
Simulations increasingly provide a number of important functions to many
scientific disciplines.5 They can be used as a phenomenological probe, suggesting new
and possibly useful directions to take in the laboratory. They can be used to test the
capability of a proposed physical model to reflect the real processes for which the model
was designed. In this paper, Monte Carlo simulations are used in order to test the
approximations used in order to obtain results from a physical model for polymer chains.
The simulations use the same model and do not require any approximations in order to
obtain results. In Chapter 4 we discuss the background of Monte Carlo simulations on
lattice models. In Chapter 5 we propose an extension to the configurational bias Monte
Carlo algorithm, an algorithm designed to obtain static properties of long chain molecules
with a high efficiency. It is interesting to note a difference between Monte Carlo and
molecular dynamics simulations. The power of molecular dynamics simulations is driven
mainly by advances in computational tools. However Monte Carlo methods also gain
performance through the development of new algorithms.
CHAPTER 2
ELEMENTS OF POLYMER STATISTICAL MECHANICS
The foundations of the statistical mechanics of polymer solutions are discussed in
the next section. In the succeeding section the modern developments in the understanding
of polymers at solid/liquid interfaces is reviewed. Particular attention is devoted to
theories of polymer adsorption.
Throughout, only the subject of linear homopolymers (i.e., non-branched and
uniform chemical composition) is treated. These objects were the first to be studied
theoretically and continue to bear relevance on more chemically complex polymers, the
study of which abounds: from dendrimers and random heteropolymers to highly
specialized helical and globular biological molecules.
Polymer Solutions
Most of the information about the conformations that polymer chains take in
solution is contained in the basic random walk problem. This problem dates to the
beginning of the century (at least). Random walks have been studied in discrete and
continuum space, however the basic results are invariant.
In this section, the discussion focuses on the equilibrium aspects of linear,
homopolymers solutions. No attempt is made to describe, for instance, the statistical
treatment of polymer dynamics, phase equilibria, or rubber elasticity, although these were
prominent and early developments in the statistical mechanics of polymers.
Historical Perspective
The first to make a complete mathematical analysis of random walks was Lord
Rayleigh6 in 1919, in the interest of describing the vibrational states of sound waves.
Einstein, Smoluchowski, and Chandrasekhar contributed to the theory of Brownian
motion through their analysis of the random walk imposed on a micron-sized particle by
thermal interactions with its solvent.' Thus the stage was set for the description of
polymer conformations in the same way. Incidentally, the phenomenon of Brownian
motion gained prominence largely because the randomness could actually be observed
with a microscope. However, it is now known that even "nanoscopic" particles such as
solute molecules diffusing through a solvent, or the particles of a gas, have the same
mechanism of motion as their microscopic cousins.' Thus it is not surprising to find that
polymer conformations possess much in common with other well-understood
phenomena. The average random walk or polymer chain has a very small size in
comparison to its fully stretched length.
In the 1930's Kuhn played a major role in the translation of the existing
knowledge about random walks into the language of polymer conformations.8 In some
sense it is an extension of the familiar ideal gas model, with the restriction that the ith
particle of the "gas" must be one unit of length apart from the (i l)th and (i + l)th
particle. Polymer models in this class are called "ideal chain" models. Such a chain is
also referred to as a "phantom chain," since the particles have no volume and can
therefore pass through each other. In the "freely jointed chain" model, the bonds are
allowed to take any angle between 0 and 360 degrees with equal probability. In 1939,
Kuhn9 pointed out that the bond angles in more realistic models for a polymer chains
(e.g., a chain with tetrahedral angles between adjacent bonds) lose all directional
correlation after only a few repeat units. Therefore, they can be "mapped" onto a freely
jointed chain by replacing the bonds with larger and fewer hypothetical bonds which
follow the true random walk. This discards information on the scale of the original bond
length and defines polymers in a very general way.
The essential feature of the ideal chain models is that they allow for the statistics
of long flexible molecules to be described in a simple way. Like the average
displacement of a Brownian particle after a certain time, the average size (e.g., distance
between ends) of an ideal chain scales as the square root of the number of steps involved:
Rend N12, (1)
where Rend is the root-mean-squared (rms) end-to-end distance. A snapshot of a 10,000
segment chain, shown in Figure 1, is reminiscent of the path of a randomly diffusing
particle.
Elastic x-ray scattering measurements on isolated (dilute) polymer chains in
solution indicated that the square root law is not completely correct.2 The average size of
the chains do scale universally regardless of the chemical units of a polymer but with a
higher exponent.
'CONF mod' using 1:3I I
CONF10000mod' using 1:3 +
200 '
150
Figure 1. The two
walk.
200 250 300 350
dimensional projection of a 10,000-segment self-avoiding
This led the way towards an important correction to the random walk. In most
solvents, the strong repulsive interactions that prevent monomer overlap have a swelling
effect on the entire conformation. This "excluded volume" effect occurs equally between
all monomers in the chain and therefore Kuhn's method of "chain mapping" cannot be
used. Flory,8 Edwardso1 and de Gennes2 each invented entirely different methods of
accounting for the effect of excluded volume. The first two are mean-field methods and
the third comes from a parallel with magnetic systems. All three of them agreed that the
correct scaling exponent is close to v = 3/5,
Rend N35. (2)
In what follows we discuss some of the simpler models mentioned above, which
are transparent enough for presentation in some detail in the space of this paper. This will
establish the elements required to discuss more advanced subjects, primarily theories of
polymer adsorption.
Ideal Chains
A random walk has a size that is proportional to N12. There are many ways to
derive the main results of the ideal chain model. The objective is to describe the average
size of a polymer coil within a dilute solution of chains. We are interested in the average
displacement between the ends of the chain with respect to all of the spatial
conformations that are possible.
Let us use the freely jointed chain model in this example. The chain consists of N
+ 1 volume-less beads joined by N bonds. The only restriction on the position, R,, of each
bead is that the bond lengths be fixed, i.e., |R, R,+1 = 1. This length, 1, is called the Kuhn
length.
The rms end-to-end distance, Rend, is a common measure of the size of a
macromolecule in solution:
Rend = (R, RN 2) (3)
where the angle brackets are used to represent the ensemble average, i.e., the average
over all configurations allowed within the model. Another equally used characteristic of
the chain size is the rms radius of gyration:
=/Z:0(R, Ro )
R2 = (N+ 1) (R -R,)2 (4)
IN R
where Rc = (, R, (N +) is the position of the chain center of "mass." The second
equality is due to a 19th century theorem of Lagrange's (see for example, reference 8,
appendix 1) which is useful in computations.
/
R,-VR,
Figure 2: Labeling of the links in a freely jointed chain.
For the moment it is worthwhile to consider a one-dimensional freely jointed
chain. In this case the displacement in the x-direction is either I or -1. For a single-link
chain, the rms end-to-end distance, Rend, is 1, since there are only two possible
conformations. Similarly, for a two-link chain, there are 22 conformations and the rms
displacement, according to Equation 3, is /[(2(22)+0(12)+2(02))/22]1/2 = 1212
In general, for an N-step walk in one dimension, we have
Rd 2K)2 N )2NN)&) + XA
end =(x- XN Ax, I Y^ (5)
where Ax, = x, x,.i. The first term in Equation 5 obviously evaluates to NI2. The second
term is zero. The average (for i # j) of any product (Ax, ) (Ax, ), no matter how distant or
near i andj are in the chain backbone, is zero because steps in the positive and negative
direction occur with equal probability. Thus, in addition to obtaining the scaling exponent
in Equation 1 we have observed that the prefactor is exactly equal to the bond length, 1.
Next we consider a two dimensional freely jointed chain. The rms end-to-end
distance in this case is the sum of the x andy components such that
Rend = (( -x ) +( N )2N = 2 ()2), (6)
Thus, the problem amounts to finding the rms projection of a single link onto the x-axis,
A link can take any angle, 0 (measured from the x-axis), with equal probability.
Therefore,
2g
cos2 Od
(( =) Y_)=2(cos22)=2 0 2 =- (7)
SdO
0
which establishes the same relation, Rend = IN12, for the two dimensional freely jointed
chain.
Consider a random walk on a simple cubic lattice. Each step of the walk is either
in the x, y or z direction. One third of the time (on average), the step is in the x-direction,
so the average squared x-component of a step is 12 / 3, much like the case of the one-
dimensional walk. In fact the presence of the lattice allows the problem to be decoupled
in terms of three independent one-dimensional random walks. This yields for the average
end-to-end distance: Rend= N(12/3+ 2/3+ 2/3)1/2 = IN12.
The above calculations are merely for illustrative purposes and demonstrate the
utility of simple random walk arguments (as well as the advantage of dimensional
arguments). The most general and concise derivation of the freely jointed chain theorem
is as follows.
R K Ar, f Ar K (Ar, )2+ (Ar, )(Ar) (8)
( =1 j=1 1=1 \I
The dot product inside the second pair of brackets in Equation 8 is equal to 12cos,,
where the angle between any two links, 0,, can take with equal probability any value
between 0 and 27r. The factor < cos 0, > characterizes the extent of correlation between
the orientation of any two links in a chain. For the freely jointed chain, this "orientation
factor" takes a value of 1 for i = j, and a value of 0 for i j. In other words
(cos,j) = 8,. This gives the result:
Rend IN12 (9)
for all dimensionalities and for any number N.
For other ideal chain models, it is not hard to imagine that the orientation factor
decays rapidly to zero11 with the separation i j Therefore, all ideal chains obey the
square root law. Examples of other models are contained in a book by Flory.12 One of
the more well-known is the wormlike chain (WLC). Instead of discrete bonds it has a
continuous backbone with a stiffness that decays exponentially along the backbone. This
stiffness is charactized by the persistence length. The only difference between ideal chain
models is a difference in their mechanism of flexibility. This model is in common usage
to describe DNA molecules.
Another model is the rotational isomeric states (RIS) model. This attempts to
incorporate realistic rotational potentials of specific polymers such as poly(methylene).
Flory has worked out the details for a number of polymers.12 The WLC, RIS and other
ideal chain models predict different prefactors in Equation 9, but all have precisely the
same power law in the limit N >> 1.
100
1 10
0 X
0 X
O X Rg good solvent
10
S 0 Rend good solvent
4 o 0.485 A Rg near theta
Sy=0.566x X Rend near theta
L 0
S--fiti
l 1 [
S0.594 fit2
Sy=0.412x
1 0.1
1 10 100 1000 10000
N
Figure 3: Scaling of rms end-to-end distance and rms radius of gyration of an isolated
chain. Data are from simulations on a simple cubic lattice in the presence of i) a good
(athermal) solvent and ii) a 0-solvent. Equations for fitted power laws are shown. The
fits use the points i) N=200, 400, 1000 and ii) N=400, 1000, 2000.
This theoretical prediction is one that is relatively easy to test. Experimentally,
light scattering can reliably measure the radius of gyration. It is possible to show11 that
the radius of gyration, Rg, of a freely jointed chain is related by a constant of order unity
to the end-to-end distance as R= 6 2Rend. Even non-ideal chains have an end-to-end
distance that is roughly two and a half of their radius of gyration. Figure 3 shows some
results that agree with the ideal chain prediction. However, the athermal self-avoiding-
walk displayed has better agreement with experimental results in general. This
discrepancy can be rectified by the incorporation of excluded volume effects to be
discussed later.
The central limit theorem of statistics: ideal chains have Gaussian statistics. The
sum of a sequence of random events will have a Gaussian distribution when the number
of events is large. That fact is usefully applied to the sum of the N link vectors of an
ideal chain, or the end-to-end distance, RN. One can start with the definition of a
Gaussian and apply the central limit theorem to the x-component of the end-to-end
distance Rx:
1 (Rx -((R ))2] (10)
P(Rx) exp 2 (10)
(T42i 2 I
where a is the statistical standard deviation, and N >> 1. The second moment of the
distribution is known from Equation 9. The x-component of the second moment is just
N12/3 which yields the condition P(R)R2d3R = NI2 /3 and resulting in o2 = Nl2. This
ultimately yields the probability distribution for a freely jointed chain
(/3 3 3R2 1 (11)
P,(R)= (27rN/2/3)-3/2 exp (11)
I-R2NI 2
The average end-to-end vector is R=0 since chains can take positive and negative
values of R, Ry, Rz with equal probability (this explains the symmetry of the distribution).
A small inaccuracy of Equation 11 is that it gives a non-zero probability for end-to-end
distances larger than the contour length IN. According to the statistics of the Gaussian
distribution, the probability of a chain size larger than N'121 (i.e., |R| /a> (3/2)1/2) is 22%.
The probability of a chain larger than the fully extended length, NI, (i.e., IR /a>
(3/2)1/2N1/2) is about 0.3% for a 6 segment chain, and 10-5% for a 17 segment chain. It
turns out that the other properties of the Gaussian and freely jointed chain distribution
also converge swiftly with increasing chain size. For example, the simulation of one-
dimensional freely jointed walks has been directly compared with the predicted Gaussian
distributions for walks ranging from 4 to 25 steps.13 To the eye, the simulation and
theoretical distributions converge for walks larger than 25.
It is important to note the range of spatial fluctuations of the chain is quite large.
The fluctuations in the distribution of end-to-end length are of order of the end-to-end
distance itself ((2 =N 2N). This "looseness" of chain conformations is an important
characteristic that makes it difficult to model real chains mathematically. The mean-field
approaches that will be discussed later neglect these fluctuations. This is thought to be an
important shortcoming of mean-field theories. However there exists an alternative. The
correct mathematical description of strongly fluctuating magnetic systems developed by
Landau and Ginzburg, was adapted by P.G. de Gennes to describe polymer
conformations.2, 14
One interesting application of Equation 11 is towards polymer elasticity. The
elastic force due to the stretching or compressing of a single chain in solution can be
calculated. If the ends are held at a fixed separation Rx, the number of conformations the
rest of the chain can take is directly proportional to P(Rx). Therefore the entropy is
related by Boltzmann's law as S/k, = InQ -3R /2N12 +const. Since ideal chains
have no enthalpic interactions, their free energy has only the entropic contribution. The
force of stretching the chain isf= -dF/dx or
3kT
f() =-a, (12)
N1
Thus it obeys Hooke's law of elasticity with stiffness that increases linearly with
temperature and decreases as the inverse of chain length. This approach has been applied
to more general elastic phenomena in rubbers, gels, and semi-dilute to concentrated
solutions. 11, 15,8, 16
The probability "cloud" of a polymer chain. Much like the Hartree-Fock theory
of atoms, one can formulate a "cloud" that represents the probability density of finding a
bead at a certain location. The cloud is radially symmetric and densest at the center much
like the s orbitals in an atom. The contributions of each bead of the chain are also
radially symmetric. In fact, the overall cloud is formed by the superposition of the cloud
of each bead. Clearly, the middle bead of the chain will contribute most highly to the
center of the overall probability cloud. Likewise, it stands to reason that the end beads
will give the most significant contributions to the outermost layers of the overall cloud. It
is helpful to have this picture in mind as we move forward.
pv(r) p(r) PN /2(r)
N/2< i
Figure 4: Probability distribution for a single macromolecule.
The probability, PN(R), that an N-step chain ends at R, is related to the
combination of such probabilities for two smaller chains to meet at that spot, PN(R)=PN-
M(R)PM(R). Then the probability cloud is the integral of this composition law over all
possible subchains:
N
p(R) = PN, (R)P (R)l f (13)
0
One immediate observation is that the average local monomer concentration is
proportional to the probability cloud. The concentration profile, or probability cloud, is
only partly understood for some problems.
Distribution functions are also universal over coarse scales. While the local
(atomic scale) distribution functions depend on the specific chemical shapes and
potentials, over coarser scales, they have universal limits. Information on the Fourier
transform of this function can be obtained from neutron scattering. As usual, predictions
for tertiary or higher order correlation functions are out of the question. Luckily, the
concentration of segments inside an ideal chain is low enough so that pairwise additivity
is valid. That leaves N2 pair correlation functions, N2/4 of which are non-identical owing
to the symmetry of linear chains. The average of these (overall pair correlation function)
can be obtained by a simple scaling argument.11, 2 The result is g(r) ~ 1/r for r << Rend ,
and g(r) ~ exp(-r/l) for r >> Rend. This derivation is based on the assumption that each
link has a Gaussian distribution of sizes and is thus on a coarser scale than even the freely
jointed chain link size. However the discussion at the end of the previous section
illustrated that even a subchain of only a few links is nearly Gaussian.
Ideal chain equation of state. If we wish to know the properties of multiple ideal
chains in a solution (i.e., dilute, semi-dilute and concentrated solutions), we can easily
derive them in an additive fashion. The entropy (and thus free energy) of the solution is
just the sum of the entropy of each of the individual chains in the solution. Just as the
segments inside an ideal chain do not interact, the segments of different chains do not
interact. For this reason (and others), it is clear that there will be much more severe
discrepancies than the difference in chain size between ideal and real polymer chains.
In practice, a polymer solution typically has a constant temperature, composition,
and volume. Thus the fundamental thermodynamic entity is the Helmholtz free energy,
F(T, V, N, n) where N represents the chain length and n represents the number of chains.
With all fundamental equations of state, any thermodynamic property is directly
available.
For ideal chains it is simple to find the difference in free energy between a
polymer solution and the equivalent systems of pure solvent and pure monomer. The
internal energy of the polymer solution is zero since all pair potentials are individually
zero at all distances. That leaves only the entropic contribution.
We will derive the equation of state for a lattice model, in a way analogous to the
Flory-Huggins method described in the next section. The number of conformations that n
chains each with N links can have is 2 = znN. Letting c = nN/V be the concentration of
monomers,
S / V =knNlnz /V kc. (14)
It is similar to the result that is obtained upon taking only the linear term in a virial
expansion. The introduction of excluded volume into the model will result in c2 and c3
terms. A fundamental equation of state can be a versatile tool in appropriating the
characteristics of homogeneous systems. In systems with a concentration gradient such
as adsorption, a "local chemical potential" based on the equation of state can be assigned
to describe the gradient.
Excluded Volume
The ideal chain models capture a great deal of information about polymer chains
in such a simple scheme. They successfully predict that a universal scaling exists for the
average chain size, suggesting that the two most fundamental properties of
macromolecules are their flexibility and length. The theoretical consideration of the
volume interactions in polymers adds another dimension to our predictive capabilities.
Is it worth the effort to account for the volume effects in a polymer coil, or does it
only have a minor effect due to the low concentration of monomers inside a Gaussian
coil? It can be shown15 that although the concentration inside a Gaussian coil is low, it is
enough to swell the coil to some degree. It should be expected that when the
concentration is forced to be higher, such as in situations of confinement or attraction
between chains and other objects, the difference between Gaussian and real chains will be
even more pronounced.
Flory's mean-field argument. Perhaps Flory was the first to understand the
significance of the role that excluded volume plays in polymer systems. He proposed the
following method15' 2 to determine the correct scaling of the chain size with molecular
weight. This method can be viewed as crude or elegant depending on one's perspective,
but his thoughts certainly did provoke a thorough analysis of the problem by other
researchers.
The method involves the minimization of the free energy of a single chain with
respect to its degree of swelling defined by --Rend2/NI2. For a swelling coefficient of 1,
the chain is ideal. To an approximation the internal energy is taken from binary
interactions only in a virial expansion as U/kT=B 2 and the concentration is
approximated as the number of segments in a sphere of radius Rend. The entropy is
assumed to be similar to that of the entropy of a stretched ideal chain. Essentially the
swelling or shrinking of the chain, as a result of binary interactions, is similar to the
stretching of an ideal chain with fixed ends which was quantified in Equation 12. This
results in
kTBN1 /2
F(a) = U(a) TS(a) = const + K 133 + -kTa2 (15)
/3 a3 2
where B is the second virial coefficient that is related to the amount of volume the
monomers exclude. It can be positive or negative depending on the strength of solvent-
monomer interactions. The result of minimization ofF with respect to a is
R9nd = /2N1/2 N3/5. (16)
It is useful to reconsider this problem in an arbitrary number of dimensions. If the
number of dimensions is d then the general result of minimizing Equation 15 is
Rnd IN3/(d+2) To understand this consider for the moment the 1-dimensional walk with
excluded volume. There are only two possible chain conformations, each with a length
IN which is exactly the scaling predicted above. This equation is correct for d < 4 in
fact. For the case d > 3 (of theoretical interest only), the monomer concentration
becomes so low that the ideal chain result N12 is recovered since there are no bead
interactions.
Edwards' mean-field theorem: the development of a self-consistent mean-field
that predicts the correct exponent for Rend in all dimensions. S. F. Edwardso1 derived the
correct scaling laws through a more quantitative method in 1965.
In the above discussion we made some mention of the similarity that polymer
chains have with the mechanism of mass transfer by diffusion. There is even a
mathematical correspondence between the differential equation for diffusion, and the
equation that describes the "diffusion" of the probability function: P(R, L). A
comparison of the mass transport and chain probability equations follows.
DV2C (17)
ot
and
dP
= 1V2P. (18)
(L
Without making too literal of an analogy, it is clear that the mathematical
solutions to the first equation can also be used for the second. Many of these have been
known for over a century from the work of Fick and Fourier.7 One can check that
Equation 11 is a solution to Equation 18. In fact, all of the results we have shown to this
point (and more) concerning random walks can be derived starting from Equation 18.
Edwards' innovation was to derive an equation like Equation 18 but with a term
that accounts for the effect of excluded volume. He solved it asymptotically for infinite
chain length to find the probability distribution:
P(R,L)= C(L)exp[ 27 (R-(R 2) 1/2 (19)
[ 20L1
where L is the contour length IN. The second moment of this distribution is
/ 2\1/2) 45 )3/ 5( j) 1 L315. (20)
3 31rl
He defined an "excluded volume," v, based on monomer pair potentials, u(r),
asv (1-e-u kT)dr .
The essential approximation in this framework is the replacement of a summation
of the "instantaneous" pair potentials at spatial position R with the average potential (or
probability) at that position. The polymer community has adopted the term used by
physical scientists (particularly in the Hartree-Fock theory of atoms): the mean-field
approximation.
It was Edwards' opinion (reference 10, pg. 614) that the mean-field
approximation was very nearly the highest level possible within the framework he was
considering.
It is a property of the self-consistent field for atoms that, although the first
approximation is straightforward though involving heavy computing, the
higher approximations are virtually impossible. The author suspects the
same situation here. It should be possible to evaluate the self-consistent
field numerically, but to improve on the self-consistent field is probably
very hard.
Indeed, since then it has required a completely separate approach, along the lines of the
Ising model in magnetics, in order to obtain better answers.
In what follows in this section we will review a simple derivation of the Edwards
equation proposed by P. G. de Gennes.2 This will prove useful since the modern theory
of polymer adsorption is an extension of Edwards' theory.
The derivation is in terms of the unnormalized probability function, or weight
function G(r ', r), for a chain on a simple cubic lattice. The weight for a chain starting at
the origin, GN(0, r), is proportional to the probability P(R, L) that the Nth step of a chain
starts at 0 and ends at r. For any conformation {ri}, the mean-field is realized as a
potential field U(r) that does not depend on {ri}; rather it depends on ({r,}). The
complete mean-field partition function is _l-Iexp[-U(r,)/kT]. Then our weighting
function is defined as
GN(r',r)=z N Z exp[-U(r,)/kT] (21)
Ir,}1, ,
where z is the lattice coordination number and the summation is only over the subset of
conformations that begin at r' and end at r, {ri}r'r. The factor zV is the total number of
conformations of N steps and can be thought of as a "partial normalization" which causes
G to have the same order of magnitude as the probability.
The weight of an N+1 segment chain is related to that of an N segment chain by
GN+(r', r) = z GN(r', r") exp[- U(r)/kT]. (22)
r
The linear expansion of the exponential term is valid since the local potential is small
compared to the thermal energy. An expansion of the weight function inside the
summation for which the derivatives vanish (due to symmetry) leads to the operating
differential equation.
dGN(r',r) U(r) a2
S G(r, r)+- V2G(r',r) (23)
dN kT 6
The partial derivative on the left follows from taking the limit as the lattice spacing
becomes small (i.e., large chain length). Equations 19 and 20 come from the asymptotic
solution of Equation 23 for infinite chain length with the appropriate boundary conditions
for an isolated chain (symmetric, etc.). A mean-field potential of U(r)~ oG(0, r)dN
based on hard sphere pair potentials was used. The approach used in the solution is the
dominance of the smallest eigenvalue for large N in an eigenfunction expansion.
Note the mathematical similarity of this to the standard Schrodinger equation
which is
-i f h 2
-ih V2+ V(r) (24)
dt 2m
The wave function t(r, t) is the square root of the probability of finding an electron at the
position r and time t. The derivative with respect to time is analogous to the derivative in
(20) with respect to N. In other words, the chain length is a time-like variable.
Integrating this over time gives the probability density, or electron cloud. Thus it is very
much like the function GN(r) and the mathematical methods can be borrowed. For a
hydrogen atom V(r) is quite simple; it is attractive and proportional to 1/r. For the rest of
the elements which have multiple electrons, a mean-field potential is used instead.
However this potential is different from U(r) in that is it based on long range, soft pair
potentials.
The Flory-Huggins Theorem A Fundamental Equation of State
A simple form for the free energy of mixing derives from the consideration of the
polymer on a lattice along with the common mean-field assumption.17'8 In the schematic
we show a square lattice, but the following applies to any lattice of coordination number
z.
F kT ln0 + (1- 0)ln(1 -0) + X( ), (25)
where 0 is the volume fraction (a3Nn/V) and X is a parameter that describes the
interactions (solvent quality). The first two terms are simply the entropy of the
configurations. The third term describes the internal energy of the interactions between
the solvent and monomers.
Notice that when N = 1, the regular solution result is recovered for the entropy.
The entropy is derived in the following way. What is the total number of
(distinguishable) ways that it is possible to fill a lattice with ns solvent molecules and n
groups of N connected monomers? In the usual way we call this number D. For a regular
solution this number is easy to find. Assuming the lattice is fully occupied by the nl and
n2 molecules of the binary mixture, there are n,,=nt+n2 total lattice spaces. The number of
ways to occupy no spaces is no!. Thus the number of distinguishable conformations
Q = n!/ln1!n,! which give the entropy, S / k = -01ln1i I2ln 2. There is no known way
to find for a set of chains on the lattice except by systematically counting all
configurations (e.g., with a computer) an impossible task for even small chains.
However, without much effort we can find a number, W, which approximates Q...
Start by adding the first monomer to the lattice. There are N1u+u possible ways to
do this for a system with n solvent molecules and n N-monomer chains. There are z-1
possible ways to add the next monomer in a way that is connected to the first. There are
approximately z-1 ways to add the third monomer and approximately z-1 ways to add the
fourth and so on. These approximations are too high especially as we add more and more
monomers to the system. A better way to estimate the number is to use an effective
coordination number z'=(z-1)(1- ). If we add the chains sequentially then the fraction of
occupied sites is changing depending on how many chains are already on the lattice 0.
In that case each chain has W,=(Ni+n )(1-0,) z',1 possibilities. Thus the product over all
1 n
chains is W =- W. The prefactor enforces counting of only the distinguishable
n!,=
chain conformations. The first two terms in (22) follow after substitution of factorials
into W, for convenience and from S = k In W.
Now we discuss the interaction energy associated with this lattice. We set the
pair potentials to be zero for all but neighbor interactions. For neighbor interactions u is
either %mm, Zsm or s. The quantity Z= X,, -(Zmm +Z,,)/2 describes the energy
difference for each monomer-solvent pair. The number of such pairs comes by
considering the number of (non-connected) neighbors of a single chain: (z-2)N+2 or
roughly zN. Then in the mean-field approximation, the number of unlike pairs is nzN(1-
0). The energy difference of mixing is then U = y(1-0). (on a per site basis) where we
have absorbed z into the X parameter.
Equation 27 was tested with osmometric measurements.2 As expected, it
performs poorly at low concentrations due to the large fluctuations present in dilute and
semi-dilute solutions. Under these conditions, W is much different than D, and U is much
different from the actual interaction energy. Thus a different form of the free energy
where the internal energy is just a virial expansion is more frequently used:
C 2 1 2 3
F/kT=-c nc+vc+wc 3, (26)
N
where v (excluded volume) and w are adjustable. This expression is in terms of
concentration c=O/v since it is meant for systems with variable excluded volume unlike
Equation 25 which implies an excluded volume of a3. A comparison of Equation 26 and
the expansion of Equation 25 yields v = a3(1-2y) and w=a3. Thus, Equation 26 has two
adjustable parameters. A solvent that has no change in internal energy upon mixing it
with a polymer has a X = 0 and is called athermal. For this case, the satisfying v = a3
results. This corresponds the typical case of a polymer in a good solvent. Further
modifications and extensions of the Flory-Huggins have also been proposed.18
Scaling Laws
A mathematical analogy between critical state magnetics (Ginzburg-Landau) and
excluded volume walk (de Gennes). It was realized that the strong fluctuations that are
typical of polymer conformations can be described mathematically in the same way that
certain strongly fluctuating magnetic systems near their critical temperature14 (Curie
point) are described. This theory was developed by P. G. de Gennes19 in 1972 and J. des
Cloiseaux20 in 1975. In the polymer analogy the Curie point is the chain length above
which universal scaling occurs. For example, in Figure 3 we have identified a Curie
point of approximately 50 lattice units. Based on this method, the most accurate
predictions for polymer chain are available. For example, the so-called Flory exponent v
(Ren,,NV) was found to be 0.592.11 This compares well with simulation results.21
Some more detailed information on the scaling of the probability distribution of
end-to-end distances can be derived. The probability of an end-to-end distance of zero is
zero and then increases rather steeply as R0.35. At large distances the probability decay is
steeper then Gaussian, more like exp(-R25). A schematic of the unnormalized probability
distributions for a Gaussian and exclude volume chains is shown in the diagram. Both of
these distributions should have a value of 1 for the fully extended chain W(Nl)=1 and 0
above NI.
While this theory predicted the same (and correct) scaling for the radius of
gyration as Edwards' self-consistent mean-field theory (SCMF), what gave it
considerable merit was its prediction of the scaling for the osmotic pressure in a semi-
dilute.
/ kT a3N9/ 4 (27)
which is asymptotic for N>>1. This is in contrast to the scaling obtained from the SCMF
(23) of r -2. Experiments confirm the predictions of (24). Perhaps a more important
result (at least for further theoretical progress) of the analogy with critical state magnetics
is the correct prediction of what is called the correlation length which is the basis of a the
accurate, versatile and simple "scaling theory."
Scaling theory. "Scaling theory" is a popular term used to describe the host of
problems that have been tackled using transparent arguments based on the results of the
magnet analogy.2' 14 The key component to this method is the assignment of an important
length scale that characterizes the range of spatial concentration fluctuations in a
polymer solution. For a dilute solution this length is simply the radius of gyration and
thus ~-N3/5. In more concentrated solutions, the overlap of chains reduces the
"correlation length." The solution is visualized as consisting of a number of independent
"blobs" of subchains. The size of the blobs is the correlation length. Each blob on
average has no interaction with other subchains and thus if the average number of
monomers in each blob nb is large enough it has the same scaling has the dilute radius of
gyration ~-nb3/5
The correlation length should shrink as the concentration increases. Its
concentration dependence by using the fact that the onset of chain overlap must occur at
* _-N4/5. Then the scaling law (~ ~m is deduced from the overlap condition
(0*)~(0*)m~Rg. Therefore m=-3/4 and since the correlation length should approximately
be the monomer size at high volume fraction:
S= 13/4. (28)
The SCMF was the first to suggest the importance of a correlation length in the
problem but there it was found ( = l1/2. The use of either of these correlation lengths in
the scaling of applied problems including polymer adsorption can produce markedly
different predictions.
Phase Diagram
This leads us to the discussion of the phase diagram of polymer solutions.4 The
length of polymer chains leads to an extra phase (semi-dilute) that does not exist in
regular solutions. The chain length and solvent quality will play the primary role in
determining the state of a polymer solution at a given concentration.
SEMIDILUTE
MARGINAL
0 DILUTE
CONCENT
1/2 IDEAL RATED
0 1
Figure 5: Sketch of the 5 regimes of a polymer solution.
The onset of the semi-dilute regime occurs when chains begin to overlap. This
"overlap concentration" c* is easily identifiable. At the overlap concentration the
solution will have the same concentration as the average concentration inside the
unperturbed polymer coils which is approximately c*=N /Rg3 N4/5. Therefore the
overlap concentration decreases rapidly with chain length and the width of the semi-
dilute phase increases. The overlap concentrations have been calculated in the chapter on
simulations for a self-avoiding walk on a simple cubic lattice. These calculations are
based on a slightly more accurate factor 47r/3 in the formula and are quite consistent with
data on the onset of shrinking of the average chain size.
One way of trying to quantify the various regimes of a polymer solution is starting
from the expanded mean-field free energy4 Equation 26. Differentiating with respect to
the number of solvent molecules and expanding the logarithmic term gives:
-p = / N+ u02/2+ ox 3 /3 (29)
The osmotic pressure is proportional to the chemical potential. Thus, the first term
represents the ideal (ideal gas) contribution, the second gives the binary excluded volume
interactions and the third term should be dominant for concentrated solutions. Five
regimes can be defined according to Equation 29. The ideal, marginal and concentrated
regimes are said to occur when the first, second and third terms respectively dominate.
The other two have already been mentioned: dilute regime where v causes a swollen
radius -N', and the semi-dilute regime, a situation of overlapping chains with a positive
v. There is a crossover between semi-dilute and marginal when the overlaps are strong
enough that the solution becomes homogeneous.
Polymer Solutions at Hard Interfaces
One of the more well established areas is the theoretical description of polymers
adsorbed to a solid surface. For the most part, solid surfaces have a net weak van der
Waals attraction with the monomer units of polymer chains. For simplicity, it is often
thought of as a square well attractive potential.
From a free energy standpoint three are three primary effects of adsorption:
reduction in the entropy, an decrease in internal energy for the segments that are
adsorbed, and an increase in internal energy as a result of monomer crowding near the
surface.
Mean-Field Theories of Adsorption on a Lattice
The mean-field lattice theory was first proposed by Roe22 and Helfand et al.23 and
reached its present stage at the hands of Levine et al.24 and Scheutjens and Fleer.25 The
main result of this work is a numerical prediction of the concentration profile of
monomer units next to a surface. It also gives the exact contribution to the concentration
profile of segments belonging to tails and loops.
The system is divided into M layers parallel to the surface with L lattice sites in
each layer. At z=0 there is an attractive (free) energy X,, and polymer concentration
decays to its bulk value as z approaches M. Edge effects will be small for large systems.
Conformations are distinguished only by the layers they pass through. Two
conformations are considered to be indistinguishable if their corresponding segments are
in the same layer, regardless of where they are located within the layer. The complete
specification of a conformation, c, would be (1,il)(2,i2)(3,i3)...((N,i) where the first
number is the segment rank within the chain and i, is the layer that segment is in.
The method is an extension of the Flory-Huggins method that was used for bulk
solutions, The canonical ensemble is chosen with the partition function Q({nc}M,L,T).
The set {nc} represents the set of conformations such that nc is the number of
conformations of the c type. Thus the total number of chains within the system is the
summation over all possible conformations: n = _nc. The volume fraction in layer i can
C
1 ,
be calculated from the distribution of conformations by =1- n,,c where nzc is the
L c
number of segments of conformation c that are in layer i.
The Helmholtz free energy is constructed from F= U kTln where U and DQ
are functions of the distribution of conformations {n,}. The change in free energy per
lattice site between a bulk Flory-Huggins solution and an adsorbed layer is actually used
AF = U- UF kTln (30)
QFH
The degeneracy DQ is constructed by adding one conformation at a time from an
arbitrary conformation distribution. The degeneracy of each successive conformation
depends only on the conformations that have already been added to the lattice. It depends
on these in an average way. In other words the probability for a new segment to be
placed in layer i depends only on the number of segments from previous chains in layer i
(Bragg-Williams approximation). Likewise the internal energy of each layer i is equal to
ZL(l- ,)( 1 + 1,+). (31)
The solution to this puzzle is to find the equilibrium distribution of
conformations. The variational principle of thermodynamics applies to the replacement
of N solvent molecules by a new chain with conformation d such that
dF = Pihandnd + pdn where upo and no are the chemical potential and number of the
solvent molecules. That results in the operating equation
(dF / d)M,L,T,nc, = chaw -Np . (32)
This yields an equation for the equilibrium number of conformations of type d,
however it is not readily solvable as such since the index d is not exactly a scalar. It is
recast in the form of the familiar probability p, that the ith segment of a chain will appear
N-1
in layer j. This yields a matrix equation of the form P,, = (w,,)) p,, for which a
numerical solution method was developed by DiMarzio and Rubin.26 The volume fraction
in each layer is then calculated by the same method as in Equation 13 and the individual
contributions are obtained analogously to Equation 41.
This theory gives comprehensive results for finite chain lengths and for
conformational features. Scaling theories, of course do not contain any information about
finite-sized chains. That means that the only other available method to obtain
information below the scaling limit is simulation. We will be making comparisons
between the two in Chapter 8.
The Ground State Dominance Solution to Continuum Mean-Field Theories
The Edwards' equation is still valid for confined systems, it requires only an
adjustment of the potential U(r). It can be written as U(z) =f(z) +J[(z)] where the first
term depends on the surface and the second term is a result of the average interaction
energies between solvent and monomer units. For a hard wall with short range (of order
of the monomer size) attraction, f,(z) can be taken as a constant in the region for 0
Elsewhere in the solution (z>l) it is zero, and for z<0 it is infinity so that nothing can
penetrate the surface. The virial expansion Equation 26 or the Flory-Huggins form
Equation 25 is often used to obtainA (z)].
The first attempts at the use of Edwards' equation to describe polymer adsorption
were by DiMarzio,27 de Gennes28 and Jones and Richmond.29 The starting point was the
adsorption of a single chain to a surface. This was a purely academic problem since it
was known that single chain adsorption effectively does not occur. That is, even in the
most dilute conditions the surface will be crowded with chains. This is in contrast to the
depletion interaction where very few chains are near the surface, and thus the single chain
analysis is often relevant. Jones and Richmond formulated the first realistic self-
consistent mean-field theory for physisorption. They limited their solution to the "ground
state" approximation, i.e., the first eigenfunction. Apparently this gives a poor
description of tails, and as such is valid for extremely long chains in which the tails are
greatly outnumbered by loops.
Scaling Description of Polymer Adsorption
In 1981 de Gennes30 used scaling arguments to characterize the concentration
profile at the interface. He proposed that there are three important regions calling them
the proximal, central and distal regions. The proximal region extends only a few segment
lengths D (which depends only on the adsorption strength) from the surface and its
characteristics are non-universal. In the central region the concentration decays with a
power law of-4/3 as
S )4/3
S a / D < z < b (33)
z +3D/4
and in the distal region it decays exponentially...
/ b= -1+ exp(-z/b) Z > b (34)
These predictions are meant to describe a semi-dilute solution in good solvent conditions.
The exponent in Equation 33 was derived from the scaling of the local correlation length
with local concentration 0 as (~ -3/4. This has been frequently referred to as the self
similarity principle.
The objective is to describe the scaling of the correlation length in the central
region. If that is possible then the concentration gradient is directly accessible. The main
length scales of the problem are the segment size a and the radius of gyration Rg. These
are relevant in the proximal and the distal regions respectively. In the central region,
neither of these are suitable for a description of the gradient of the correlation length in
the central region. The correlation length must have an increasing scaling law with the
distance from the surface z since for adsorption, the concentration decreases smoothly
from the surface to the bulk. Near the beginning of the central region the correlation
length will be of order of the monomer size and at the end of the central region it will
have reached the bulk correlation length or in the case of the a dilute solution, the bulk
radius of gyration. From a dimensional argument alone, the correlation length must be
linear in z since no other appropriate lengths exist.
z (35)
To summarize then the resulting predictions for the concentration profile:
DeGennes ~ -3/4 0 z -4/3
Edwards 0-1/2 o~ Z-2 (36)
0 solvent 0-1 -- 0 ~ z-1
These have not been tested conclusively. The second one listed should be appropriate for
the marginal regime. The conditions inside the central region usually fall within the good
solvent / semi-dilute regime, therefore the first part of Equation 36 is the most useful.
Square-Gradient Mean-Field Free Energy Functional
In the same paper de Gennes developed a mean-field approach that is an
extension of the Cahn-Hilliard theory for the free energy of a region of regular solution
near an interface. This theory proposes that concentration gradients in a solution add a
term proportional to the square of the gradient in the free energy expression. Thus his
expression for the free energy of a polymer solution interface, for a 0-solvent is
Y Yo= r, Y + d-- ++ dz (37)
where V is a complex number that determines the concentration profile by 0 = I V 12. The
first term is the contribution from the surface adsorption, the second is related to the
entropy of a concentration gradient (as mentioned above), and the third accounts for the
interaction energies between monomers (see Equation 26). Thus for any volume fraction
profile, (z), the difference in free energy between the interface and the bulk is fully
determined. The (z) that produces the lowest possible free energy difference
corresponds to equilibrium. It is interesting that the minimization produces the same
differential equations that would result from the Edwards equation. An approximate
analytical solution was found for the central region in a 0-solvent38 as
0 a (38)
2co(z+ D)
therefore the expectation is that a broader layer is produced in 0 -solvents. Inherent in
this solution is ground state dominance. Finally, in this comprehensive work, a slightly
different form for the free energy in Equation 37 is suggested that can produce the scaling
results Equation 33 through a minimization procedure identical to the one used in the
mean-field approach.
Beyond the Ground State
Ploehn and co-workers31 extended the mean-field theory of Jones and
Richmond.29 They attempted to accurately model the boundary proximall region) with
the idea this might have an impact on the overall structure of the interface. A second
"boundary Edwards equation" is matched with the regular Edwards equation. This
system was numerically solved for various chain lengths and solvent qualities. The
choice of the potential derived indirectly from Flory-Huggins Equation 25:
U(z)/kT= -In (-) -2X( Ob) (39)
(1 -O b)
Although solutions for various chain lengths were obtained, only those for the highest
chain lengths should be correct due to the ground state dominance approximation.
Several developments have come recently from Semenov, Joanny and co-
workers.32, 33, 34 Using the same square gradient approach as de Gennes but with an
interaction energy in (23) dominated by v, they find
O(z) = Ob th2(z / 2 + const) -o > 2/(z+ D)2. (40)
This is supposed to be appropriate for the so-called marginal solvent conditions.
Separately they start from the Edwards equation and attempt to find a solution
that takes into account the effect of tails. Two eigenfunctions are used to describe bound
and free states. The partition function G(n,z) is divided into adsorbed and free parts Ga
and Gf. The contribution that loops and tails make to the concentration can then be
calculated by
01(z) = dnG, (n,z)Ga(N -n,z)
(41)
(z) = 2b f dnGf(n,z)G (N-n,z)
N
Additionally they proposed from scaling arguments that the central region can be split
even further into two regions: call them central a and central b. Central a is dominated
by loops and central b by tails. The crossover in the concentration profiles z* scales with
the chain length as N1/3
Juvekar et al.35 followed up on the work of Ploehn et al.31 recently by extending
the numerical results to a system with variable solvent quality. Again they take an extra
boundary Edwards equation and even a boundary solvent quality in their solution. They
suggest that a simplification used previously in the mathematics was unwarranted. There
have been indications that the "effective" solvent quality can change with chain length
and concentration. They take this effective X to be of the form:
= xZoO( +c+ c202)(+ C3/ N) (42)
They suggest that the effect of this improves agreement with experimental data for bound
fraction and adsorbed amount. Also they can calculate the tail and loop profiles.
CHAPTER 3
EXPERIMENTAL CHARACTERIZATION OF POLYMER INTERFACES
Measurement of Polymer Structure at Surfaces
Knowledge of the structure formed by adsorbing polymers at interfaces is
essential to the ultimate understanding and control such layers for engineering and other
purposes. It is a demanding experimental task to measure, for example, the structure of a
layer of physisorbed homopolymers at equilibrium with a bulk solution. This requires
carrying the experiment out in a liquid environment and the use of non-invasive
techniques. These caveats limit the options considerably. For example, even x-rays exert
some physical influence on polymer layers in some reflectivity or scattering experiments.
Physical probes, even in the rare cases when they are as "gentle" as the AFM tip in
tapping-mode or hydrodynamic probes, cannot image the internal regions of a polymer
layer.
One of the few approaches that meets the requirements of such proposed polymer
measurements is the utilization of neutrons. These have a small enough wavelength to
probe the atomic scale and are completely non-destructive owing to their infrequency of
collision with sample atomic nuclei.
In the following section we discuss some of the overall, or global properties that
can be measured without resort to a neutron source. Subsequently we discuss the detailed
structural information that only neutron scattering and reflectivity can probe, and their
current limitations in sensitivity for adsorbed layers. An extensive review of all of these
techniques is contained in the Chapter 3 of the book by Fleer et al.4
Global Properties of Physically Adsorbed Layers
There are several static properties of polymer chains at an interface that can be
measured experimentally. They consist mainly of three types: measurements on the total
adsorbed mass, layer thickness, and the fraction of monomers in a chain that are directly
bound to the surface.
The adsorbed mass, F, represents the total number of chains in contact with the
(i.e., adsorbed) surface per unit area, multiplied by their molecular weight. Thus its units
are in terms of monomers per unit area. An equivalent quantity can be obtained from
theoretical or simulation data for the adsorbed-chain volume-fraction profile, 0a(z), as
F = J (z)dz= (z). (43)
0 z=1
The summation is written for use with theoretical models based on a lattice, and
consequently with a discrete concentration profile.
The adsorbed amount can be measured by flushing the solvent and excess
polymer from a system in which the chains have had sufficient time to adsorb. The
adsorbed amount is the difference between the weights before and after adsorption.
Provided desorption kinetics are slower than the flushing process, this is usually a reliable
method.4
Layer thickness is an ambiguous term as the adsorbed layer is diffuse and does
not have an abrupt "edge." Layer thicknesses can be defined in several ways. The mean
distance of polymer segments (belonging to adsorbed chains) from the surface, or for that
matter, any of the statistical moments. The most frequently used actually is the second
moment, or rms layer thickness:
SJz 20(z)dz
8.0 (44)
Again a discreet volume fraction profile obtained from simulations would require a
summation instead of the integral.
While there is no technique that actually measures this hypothetical quantity there
are some techniques which are thought to contain related information. In a dilute
dispersion of spherical particles, neutron scattering is said to measure the difference
between the squares of the second and first moments.4
One of the most common and accurate ways to measure layer thickness is by light
ellipsometry (similar to reflectometry). The sample is a macroscopically flat surface
covered with adsorbed polymer inside a fluid cell. Light is reflected off of the surface at
a low angle and the intensity distribution of the reflections give information about the
refractive index of the polymer layer as a whole. There have been several attempts to
link the quantity that is measured by ellipsometry to the theoretical volume-fraction
profile. The definition preferred by Fleer et al.4 and which we will use in succeeding
chapters is
6ell f0d2 (45)
The most intuitive quantification of thickness is through hydrodynamic methods
because these quantify the degree to which a fluid flow is influenced by the presence of
the layer (as compared with flow past the bare surface). Fluid is driven past a polymer
layer adsorbed to the surface at z = 0. The extent to which the flow is reduced is
interpreted as equivalent to the linear flow past a bare surface at z -= hyd. Typically the
polymer is adsorbed on the inside of a capillary tube. The Debye-Brinkman
hydrodynamic model for flow past a porous layer can establish the relationship of the
hydrodynamic thickness with the volume fraction profile.4 The velocity profile in the
direction normal to the surface v(z) is described by the following equation of motion:
d2v _
=v (46)
dz c(1- (a)
where the drag term on the right depends on the relative porosity of the layer (the
permeability constant c is usually taken as 1.0. The linear part of solution for v(z) is
extrapolated to zero and the z to which this corresponds is taken as the thickness.
The bound fraction p is the fraction of the segments belonging to adsorbed chains
that are in contact with the surface. Several methods exist to observe this, for example
infra-red spectroscopy can distinguish certain types of contacts with the surface such as
hydrogen bonding. Another method that can distinguish between attached segments and
non-attached segments in nuclear magnetic resonance. A magnetic pulse is applied to the
system and the relaxation of the polarization of the atomic nuclei gives information about
the thermal mobility of these atoms. Presumably the thermal motion of attached
segments is greatly reduced. The bound fraction is easy to obtain from lattice
calculations since the number of segments in layer 1 is just p = (1) / F.
Small Angle Neutron Scattering
There exist few experimental methods with sufficient sensitivity to probe the
polymer volume-fraction profile. For reasons stated above the most progress has been
made by small angle neutron scattering (SANS) and neutron reflectivity. An excellent, if
not amusing, account of neutron scatting that assumes little foreknowledge of optics is
available in the review by R. Pynn.36 Reviews of neutron reflectivity experiments on
polymers have been written by T. P. Russell.37' 38 A non-technical overview of these
techniques follows.
The length scale of neutrons and x-rays, of the order of 10-7 to 10-10 m, is
appropriate for probing the microstructure of many complex fluids. Neutrons follow the
same quantum mechanical principles as do photons and thus neutron scattering has much
in common with conventional light scattering. In fact the refractive index has a parallel in
neutron techniques which depends on the makeup of a sample's atomic nuclei the
scatterers. A measure of the strength of neutron-nucleus interaction is the scattering
length, b. For example the coherent scattering length for hydrogen and its isotope,
deuterium, are -3.74 and 6.67 fm respectively. Actually this large difference turns out to
be very important in contrast matching. By adjusting the level of isotopes in a material,
its refractive index can be controlled. Thus the feature of interest can be labeled by
adjusting the contrast of the other components in the system so that they are "invisible."
A SANS experiment on polymer layers proceeds as follows. A dilute dispersion
of colloidal particles bearing polymer layers is prepared, using deuterium isotopes to
match the refractive index of the solvent to that of the colloids. The key parameter is the
scattering vector, Q, which explores the range over which structure can be explored. The
scattering vector is varied by adjustment of the wavelength, A, and the detector angle, 0,
defined by
47y 8A
S= sin (47)
>A 2
where the neutron refractive index, p, is a derivative of the scattering lengths in the
materials.
The output of the instrument is the scattering intensity as a function of the
scattering vector, I(Q). Models have been developed to interpret the structural
information in I(Q) which describe the scattering due a thin spherical shell (i.e., the
polymer layer). Thus the scattering intensity is mathematically related to the volume
fraction profile through an integral. In one such model the limit of low Q yields4
const !
I(Q)- os (z)ezdz (48)
Q 0
At this point a difficulty in the analysis becomes apparent. The inversion of the integral in
Equation 48 to obtain (z) is not trivial. There have generally been two methods used to
proceed. The first method uses assumptions about the functional nature of (z) and a self-
consistent fitting procedure.39 The second method does not invert the integral in Equation
48 at all but rather analyzes the data directly for certain features that would be consistent
with predicted layer structure.40
Neutron Reflectivity
Neutron reflectivity is similar to neutron scattering. The major difference is that
instead of a colloidal dispersion, the experiment is carried out on a macroscopically flat
substrate bearing the polymer layer. The effectiveness of the reflectivity technique has
emerged in the last 5-10 years. Its major strength is its high depth resolution to
approximately 0.5 nm.38
The analysis is similar (broadly speaking) to that described above for neutron
scattering. Reflections at low angles off of the polymer layer and substrate can be
analyzed in terms of a sequence of Fresnel films of decreasing refractive index for the
case of homopolymer adsorption. Reflectivity has found specific merit in the analysis of
diblock copolymer interfaces.38
Scattering and Reflectivity Measurements on Adsorbed Polymers
Relatively few SANS41' 42, 43 and neutron reflectivity44' 45 studies have been
performed on physisorbed homopolymers layers in a good solvent. The SANS
experiments of Auvray and Cotton42 (polydimethylsiloxane (PDMS) with MW=270,000
adsorbed on silica) were consistent with the scaling power law prediction for the
intermediate interface (-4/3, as opposed to -2, the mean-field prediction). However,
neutron reflectivity measurements on shorter deuterated polystyrene (PS) (MW=50,000)
adsorbed on mica were consistent with a simple exponential decay of polymer
concentration.45 However the sensitivity of these techniques is such that volume-fractions
below 1% cannot be observed,39 which constitutes a substantial part of the profile from
physically adsorbing chains (unlike grafted chains).
It is notable that several attempts have been made to observe simultaneously46' 47
or in parallel48 the structural aspects and the interaction forces between polymer layers.
This is certainly a direction that will lead to a richer understanding of the cause and effect
of the response of polymer layers. A new surface force apparatus which uses neutron
reflectivity in situ is described47 and force and volume-fraction profiles for polystyrene
adsorbed on quartz in cyclohexane were obtained. In a similar experiment using block
copolymers the volume fraction profile was observed to increase due to the compression
and a collapse in the layers was witnessed when the solvent was changed from good to
poor.
Measurement of Polymer Induced Forces
The primary tool that has been used to measure forces between surfaces bearing
polymers is the surface force apparatus (SFA) developed largely by J. Israelachvili.49' 50
The SFA remains the only apparatus capable of measuring directly the forces induces by
polymer chains between semi-infinite solid surfaces. For this reason the bulk of this
section will discuss the apparatus itself, and measurement that have been made on
homopolymers physisorbed from a good solvent. Much progress in this area was made in
the 1980's, however in last decade SFA studies have focused on more complicated
systems.
Another tool that has great potential towards isolating the forces induced by
polymer chains is the atomic force microscope (AFM).51' 52 While the AFM uses an
inherently microscopic, or local, probe, continuing growth in this area has led to advances
such as the use of spherical 10 jam particles as probes.53 One advantage of AFM is the
variety of materials that can be used for the two solid surfaces in contact. In the future,
advances in AFM will complement the understanding gained with the SFA.
Other techniques for measuring surface forces have been reviewed by Claesson et
al.52 but, in general, have limited applicability to polymer systems. Total internal
reflection microscopy on a colloidal probe controlled by a laser has been used to measure
weak forces of order 10-14 N. Forces can also be inferred indirectly from phenomena such
as the swelling of clay minerals.54 Naturally prior experience with the phenomenon as a
"calibration" of the technique is prerequisite to these approaches.
Surface Force Apparatus
The surface force apparatus is an instrument that can sensitively (10 nN, 0.1 nm)
measure the force and distance between two curved mica surfaces, often enclosed in a
liquid cell. This is accomplished with a combination of spring deflection and multiple
beam interferometry. A cursory description of the apparatus is given here as a reminder
of its basic function and capability. For an extensive discussion of the development and
details of the SFA the subject is well reviewed.52' 55 3, 56
The major success of the SFA was the confirmation DLVO theory. The use of
molecularly smooth mica surfaces was instrumental in factoring out the substantial
effects of surface roughness in van der Waals forces. Mica however is also important in
the method for determining the distance between the two surfaces because of its optical
transparency. White light passed through the specially coated surfaces yields the distance
through interpretation of the interference "fringe" patterns associated with the reflection
of light off of each surface.
The fine separation between the surfaces is controlled by a piezoelectric ceramic.
Each volt of electricity applied induces a 1 nm expansion of the tube-shaped piezo.
Ultimately the force is obtained by knowledge of the spring constant and the difference
between the actual and applied displacement. Clearly then, reliable measurement of the
separation is doubly important.
The mica sheets are mounted on crossed cylinders in order to avoid the difficulty
of aligning the surfaces in parallel. Force profiles are always normalized by the average
radius of curvature of the cylinders so that they are independent of geometry.
I Piezo tube
__ Leaf spring
Light
Figure 6. Elements of the surface force apparatus. Light passes through the lower mica-
coated hemi-cylinder which is mounted on a leaf spring in order to measure the force.
The silvered surfaces of the hemi-cylinders produce an interference pattern that is
collected above. The upper cylinder is mounted onto a calibrated piezo-electric ceramic
which enables fine adjustment of the distance between surfaces.
Usually a polymer adsorption experiment will proceed as follows. A polymer
solution is introduced into the liquid cell the system is allowed to "incubate" for a period
of approximately 12 hours. When the incubation, or adsorption, is complete, the liquid
cell is washed out with pure solvent and the compression / decompression cycles begin.
Typically the rate of compression is 5 to 60 minutes. This is supposed to be a quasi-
static, or equilibrium, process.
SFA Measurements in Good Solvent
A number of measurements have been made on physisorbed homopolymers in
good solvents.57' 58, 59, 60, 61, 62, 63, 64 These experiments as well as measurements conducted
in poor solvents have been reviewed comprehensively by Patel and Tirrell,65 and also by
Luckham.66 The good solvent systems that have been studied are listed in Table I.
Table I. Experimental characteristics.
Polymer Solvent MW Reference
PEO Aqueous 40,000 57, 58, 59,
Electrolyte 160,000 60, 63
1,200,000
PEO Toluene 40,000 61, 62
160,000
PS Toluene 65
Of those that are in regular usage, there are not many polymers that adsorb onto
mica surfaces from a good solvent. Polyethyleneoxide (PEO) is the only one that adsorbs
strongly enough to make consistent measurements.61 65
In all experiments where the polymer was given sufficient time to adsorb, the
forces observed were monotonically repulsive, growing in strength with decreasing
separation. By contrast, all observations in poor solvents shown an strong attractive well
with its minimum at a separation of roughly 1 Rg.67' 68 The range of interactions in good
solvents is less exact, but generally corresponds to 5 10 Rg.59' 61 For poor solvents
interactions occur only for only 2-3 Rg, indicating the swollen layers found in good
solvents. Toluene is a better solvent than water for PEO and this results in an onset of
repulsive forces several Rg before the aqueous system.
Clearly if polymers are bridging between the surfaces, this will lead to a (probably
substantial) attractive contribution to the force. There has been some uncertainty about
the role of bridging. The simulation study by Jimenez et al.69 indicates that bridging
occurs in all systems, and the number of bridges does not vary much, except with
separation. The osmotic pressure is therefore the root of the repulsive forces observed in
good solvent systems. The monomer-solvent potentials in good solvents are such that
monomers are effectively repulsed from one another due to the entropic, and in very good
solvent, the enthalpic benefit of solvent contact. Therefore the osmotic forces which are
attractive in poor solvents, are strongly repulsive in good solvents to the extent that they
can overcome contending bridging forces.
The largest molecular weight chains used, 1,200,000 Mw PEO in water could be
sufficient70 to test asymptotic theoretical predictions provided polydispersity
(Mw/Mn=1.2) or other effects do not play a substantial role. The density functional
approach used by P.-G. de Gennes30 yields two power laws for the force profile at far and
near separations. The data are consistent with the inner power law (F D-2.25) but an
outer power law (F D-3 is difficult to discern. However the authors preferred to
interpret the results using an analogy between the tail conformations of a strongly
adsorbed layer and a grafted layer. This was motivated by the fact that the force profiles
a common qualitative feature with grafted polymer measurements when shown on a
semi-log plot. On approach of the surfaces the force rises sharply, flattens briefly and
rises steeply again. However, the role of tails should dissipate for exceptionally large
molecular weights and adsorption from a dilute solution. Furthermore, the hysteric
effects noted on decompression may be indicative of weak adsorption.
Klein and Luckham have studied the effect of undersaturation on the forces in a
PEO/aq. System.59' 63 An attractive well of depth -50 [IN / m is observed for such low
surface coverages. This is small in comparison to poor solvents (by a factor of 20) but is
probably enough to explain the flocculation of colloids. This is explained by the relative
reduction in osmotic forces since fewer monomers are present between the surfaces. The
method employed in generating undersaturated layers is as follows. The surfaces are
brought close together (several microns) before the polymer solution is introduced. This
creates a diffusion limited adsorption process. Thus, long times are required for the
layers to become fully saturated. Compression cycles are performed at several intervals
during the slow (i.e., 48-hour) adsorption process, each corresponding to a different level
of undersaturation. With increasing surface coverage, as the minimum becomes
shallower, it also shifts towards larger separations and ultimately disappears well before
full coverage is reached. This is consistent qualitatively with theoretical models
including mean-field,71' 72 scaling73 and simulation69 results, although between these there
are quantitative differences for the range of undersaturation that produces attractive
forces.
Many researchers have noticed considerable hysteresis in their measurements in
good solvents. Unresolved differences in hysteresis using the same polymer / solvent
have even been reported.59' 62 In general good solvent measurements have been less
reproducible and well behaved than poor solvent or grafted polymer measurements.65
There are several conceivable reasons for differences between decompression and
compression force profiles. Perhaps constrained equilibrium (fixed number of polymer
chains) is not plausible for slow compressions, since they could diffuse from the area of
contact. This might effect would probably be strongest for large surface coverages74 and
weakly adsorbing surfaces. Another possibility is that the compression, in compacting
the layers, forces many segments to adsorb and the desorption / relaxation of these
compact conformations is slower than the decompression. Weaker forces on
decompression have frequently been observed.63 In one experiment,59 after a rapid
compression, a constant force was applied and the surfaces decompressed at their own
speed. This produced the same force profile observed from a slow (equilibrium)
compression. Little progress has been made toward robust explanations and predictions
of non-equilibrium effects in these measurements. This is an area where dynamic
computer simulations could make a substantial contribution.
Atomic Force Microscopy
The first AFM was developed by Binnig, Quate and Gerber.75 In addition to the
impressive ability to measure topographical detail of surfaces to the order of tenths of an
angstrom (e.g., the well known imaging of the hexagonal structure of carbon-carbon
bonds in graphite76), the apparatus has been outfitted to measure forces parallel and
normal to a surface. Thus it has contributed to the growing field of tribology (atomic
scale friction) and measurement of surface interactions. Variations on the basic
instrument abound and its uses have grown to include the semi-conductor industry and
biological applications, in addition to atomic and electronic structure and interactions.
The basic instrument consists of a tip or probe attached to a cantilever-spring and
a sensor to detect the deflection of the cantilever. Supplemental components include
piezo-electric system for moving the tip laterally and normally, and a feedback system to
control the applied force, or mode of the tip. Generally there two classes of operating
modes: those that image a surface and those that measure forces between objects.
A constant force "contact mode" is used typically for imaging of hard surfaces.
The high resolution (subangstrom) of this mode is achieved by using a small loading
force between 10-7 and 1011 N.51 The "tapping mode" is typically used for imaging soft
surfaces. measurement of long range forces This uses the method of "resonance shift
detection." The cantilever vibration is driven at its resonance frequency. The vibration is
substantially dampened (or amplified) in the presence of long range repulsive (or
attractive) forces. This enables the instrument to achieve even smaller loading forces of
order 10-13 N (51) which do not disrupt the sample. Finally, and emerging method for
imaging is achieved via magnetization of the probe. An image in generated from the
response of the tip to the magnetic field of the sample.
Laser a)a
deflection
sensor b)
Cantilever
....... Piezoelectric
controlled
stage
Figure 7. AFM experiments on forces induced by polymer chains, a) An example setup
of the AFM for measuring forces between polymer layers adsorbed to colloidal particles.
b) A typical experiment on a single chain with its ends grafted to the AFM tip and the
substrate.
The above modes are used in a lateral scan of the sample. In order to measure a
the normal force profile between surfaces the technique becomes much simpler. In this
method, the voltage to the piezo is ramped (linearly) to drive the surfaces together. Often
a colloidal probe is attached to the cantilever in order to avoid geometrical complexities
in the analysis. This method was first used to measure electrostatic interactions between
colloids.53
AFM Measurements on Elasticity of Single Chains and Forces Between Polymer Layers
The strength of steric and other polymer induced forces (10-9 1011 N) has been
measured between colloids bearing adsorbed polyelectrolyte layers77 and between
colloids bearing adsorbed, aqueous PEO (MW=56, 205 and 685- 103).78, 79 In the most
recent study of Braithwaite and Luckham,79 the force profile was obtained when one or
both surfaces bear a polymer layer at a range of coverages (however both surfaces were
attractive to PEO). Mostly repulsive interactions were observed as the rate of approach of
the surfaces was relatively rapid. The range of interaction between saturated surfaces
scaled with molecular weight along the same lines as predicted theoretically
(experimental onset of interaction: N0.4, theoretical layer thickness: N0.6). No
comparison with the theoretical force profile power law30 was attempted.
Single-molecule polymer experiments are the most exciting new applications of
the AFM. These experiments started with biological molecules such as DNA and Titin
since are exceptionally large and easy to manipulate on an individual basis. In one such
experiment, a DNA molecule was covalently bonded by either end to the glass surface
and the colloidal probe.80 The adhesive force during retraction increased until a rupture
occurred between 0.8 and 1.5 nN. In a other experiments81 researchers actually observed
stretching of a DNA molecule to nearly twice its contour length indicating a rupture of
the its basic helical structure. Dextran also underwent a distinct conformational transition
during stretching at similar forces.82 Titin, a molecule found in muscle tissue, had a
sawtooth shape in its extensional force profile attributed to the unfolding of globular
units.83
Very recently, some success in similar experiments with synthetic molecules was
announced. Most notably, poly(methacrylic acid) was grafted to a gold surface treated
with surfactants in order to keep the grafting density low.81 The single molecules were
stretched by virtue of their physical attraction to the AFM tip and the elastic force profile
was fitted to the freely jointed and wormlike chain models of entropic elasticity. The
persistence length was found to be of the order of the monomer size. Also similar
experiments with polyelectrolyte84 and block copolymer85 chains have been performed.
One experiment imaged isolated adsorbed polyelectrolytes using "tapping-mode" to
obtain the lateral and normal dimensions of the adsorbed chain for different molecular
weights.86 In all, more than 30 experiments (for a summary see reference 81) have been
reported in the last 5 years that analyze the stretching of a single chain, and a new journal
has appeared entitled "Single Molecules" which takes substantial contribution from AFM
research.
CHAPTER 4
MONTE CARLO SIMULATION OF POLYMER CHAINS ON A LATTICE
Here we present an overview of lattice Monte Carlo techniques that have been
used for the simulation of polymers. Our main focus will be on the challenges and
practical considerations in our lattice simulations of polymers. A comprehensive review
of all early methods of lattice simulations of polymers was done by Binder and Kremer.87
Review has also been made of the more recent advances in Monte Carlo (and molecular
dynamics) simulation of polymers.88 5
Why are lattices still used to model polymers even now that we have such
comparatively powerful computers at our disposal? One reason is that the number of
important configurations that a polymer can take, even on a simple lattice is
overwhelming even for a modern computer. To attempt to incorporate more detail into
the model requires a sacrifice in the range of systems one can study. The classical
example of this is the scaling in dilute solution of the radius of gyration with chain
length. While atomic details of chain conformations depend on the model used, the
overall chain size is one of many universal features of polymer conformations. It would
be wasted effort to attempt to obtain the Flory exponent, v = 3/5, from a continuum
model. For the long chain lengths necessary to obtain the exponent accurately, the effect
of the lattice vanishes. Such long chain lengths are not well suited to "weighty"
continuum approaches. Therefore the most exact Monte Carlo calculations of universal
exponents have come from lattice studies.21 There are many other problems in which
coarse-grained information is desired and the advantage of these methods is evident.
Enumeration, an Exact Theoretical Treatment for Model Oligomers
In many ways, Monte Carlo simulation is a substitute for the full enumeration of
all possible configurations of a system. On a simple cubic lattice, a 100-step random walk
can trace 6100 1078 conformations. While the fraction of these that are self avoiding is
extremely small (10-26 %), the number of self avoiding walks is nevertheless staggering.
For the same size chain a self-avoiding walk (SAW) has ~1050 conformations. Current
high-speed computers can only generate of order 109 conformations,4 thus a modern
computer cannot even begin to enumerate the statistics for polymer problems of modern
interest. Enumerations are not simulations, they are exact statistical mechanics
computations (for lattice models).
Enumerations seek to generate each term in the summations representative of
thermodynamic quantities. For such a property, A, in the canonical ensemble, its standard
configurational average is
A(rN)exp[-P(f(rN)]
(A)NV = rN (49)
rN
where the summations range over the entire set of states possible for N-components on a
lattice of size V. Actually the denominator in Equation 49 is often enough to calculate
properties of interest.
Consider the exact enumeration of a SAW of only two steps. On the simple cubic
lattice with coordination number z = 6 and step length a, there are
S= z(z -1)= 65 = 30 self-avoiding walks that are possible once the first bead has been
placed. One-fifth of these are straight conformations with an end-to-end distance of 2a,
while the other four-fifths are bent conformations with an end-to-end distance of 212a.
Therefore the average end-to-end distance is -2a + 2a = 1.55a .
The number of self-avoiding walks grows quite rapidly. For long chains this
number has the following scaling.2
Q,, = cz"N (50)
where y (=1/6 in 3 dimensions) is a universal exponent called the susceptibility, c is a
constant and z ( < z )is the effective coordination number roughly equal to 4.5 for the
simple cubic lattice. Enumeration studies are usually limited to chains with N 20.87
As the number of steps, N, increases, the fraction of all random walks that are
self-avoiding is vanishingly small. This is evident from the limit of the ratio
Q =A c /N16 (51)
Q RW Z
which approaches zero for large N.
Monte Carlo Methods
Static Monte Carlo
Monte Carlo integration. The first Monte Carlo methods were used for relatively
simple numerical integration, or convolution, as it is called. Instead of using techniques
like Euler's method to evaluate the integrand at a sequence of points within the range of
integration, in Monte Carlo integration the integrand is evaluated at a large number of
random points. This proves to be useful for quickly varying functions.
The advent of statistical mechanics made it necessary to improve the methods of
multi-dimensional integration. To compute integrals or summations like Equation 49, a
microstate, r", must be generated at random. The attempt must be abandoned if, in the
SAW example, there is a multiply occupied lattice site. Even for a single chain, though
not quite as severe as Equation 51 (because only non-reversal random walks (NRRW) are
considered), the number of successful attempts is discouragingly low. Thus little or no
progress was possible in more concentrated systems.
Importance sampling. Improvements were made on the Monte Carlo method that
sample the integrand in a range biased towards the most important contributions to the
integral. "Importance sampling" must be corrected with the appropriate weighting factor
for the non-uniform distribution from which it samples. The improvement of the
sampling distribution is the focus of ongoing development in Monte Carlo algorithms.
Of many proposed static Monte Carlo methods for polymer chains (e.g., simple
sampling, dimerization, enrichment), the most successful is the biased sampling proposed
by Rosenbluth and Rosenbluth89 in 1955. In this method, a chain is generated one step at
a time. The chain avoids overlaps during its growth by selecting from only the
unoccupied neighboring sites. They were able to induce, from a simple example, the form
of the bias that this selection procedure imposes.
Consider the two-dimensional construction of a 4-step walk. There are exactly
100 possible 4-step SAW's starting from the origin. Consider the following two defined
by coordinates (0,0) (1,0) (1,1) (0,1) (0,2) and (0,0) (1,0) (1,1) (2,1) (2,2), shown in
Figure 8. These should occur with equal probability since we are not assigning any
interaction energies. However in this scheme, move A occurs with a probability
(1/4)(1/3)2(1/2) while move B occurs with probability (1/4)(1/3)3. The properties of A
need to be weighted by the factor (2/3) in order to correct for the bias.
A B
Figure 8. Two of the 100 conformations available to a 4-step self-avoiding walk and a
two dimensional square lattice.
To generalize, the weighting function, WN, required to correct any given N-step
walk is defined by the recursive relation (W1 = 1):
W+1 = (n /5) W,. (52)
where n is the number of empty neighboring sites at the mth stage of the walk. This
weighting function is used for thermodynamic averages, instead of Equation 49 we have
(for the case of zero energy of interaction)
C. W,(,-)A(r-")
(A)N =im (53)
where the summations range over the sample states generated (and therefore it is an
approximation). This algorithm generally works well for single chains smaller than 100
segments. Its power is considerably extended when substituted into a dynamic Monte
Carlo scheme
Dynamic Monte Carlo
The dynamic Metropolis Monte Carlo method was developed by Metropolis et
al.90 in 1953 and has been applied to everything from sub-atomic physics to colloidal
suspensions. In one sense the method is a leap forward in selection of a biased sample.
This is because it samples from a Markov chain of states. In other words, each trial state
is generated from the previous state of the system. This allows incremental changes to
occur in large systems rather than abandoning most of the trial attempts as in the static
Monte Carlo methods.
The method allows any set of "dynamic" moves to be used, provided that they
follow two important criteria: ergodicity and semi-detailed balance. A move is ergodic as
long as every configuration of the system can, in principle, be reached after a finite
number of moves. Non-ergodic algorithms can be alternated with ergodic ones which
makes the set of moves as a whole ergodic.
The semi-detailed balance condition requires that the frequency of moves any
state, i, should be equal to the frequency of moves out of that state into other states.
Normally this condition is hard to formalize and the condition of detailed balance is used
instead (a more restrictive criterion than semi-detailed balance). According to this, the
frequency of forward and reverse exchange between a pair of states, i and j, should be
identical. The condition is then subject to the transition probability, H,,,between two
states by a given move. Let P, be the Boltzmann weight, exp(-1A(r,)), of state i. The
detailed balance and ergodic conditions ensure that the number of configurations of i that
are generated will be proportional to B,. Thus the detailed balance condition is
Pn,, = Pn, ,. (54)
When a trial, j, move is fully generated, it can be subjected to an acceptance condition. It
will be accepted with a probability A,. When the move probabilities are symmetric
between all pairs the acceptance test is the only constituent of lH,. In this case the best
solution to Equation 54 is to use the following acceptance rule.
A, = min [ 1, B / B, ]
(55)
Therefore a Metropolis scheme is first proposed and the appropriate acceptance rule is
derived from Equation 55. This forces the sampling distribution to be the Boltzmann
distribution. The consequence of this is that the Boltzmann weight drops out of the
thermodynamic averages of Equation 49 yielding instead
(A) = lim (56)
One of the earliest algorithms designed for lattice polymers uses local rotations of
one or two segments along the backbone of the chain called kink-jump91 and crankshaft
moves."7 This is currently an excellent method to use for melt-like concentrations. Also
its correspondence with Rouse dynamics makes it attractive for kinetics studies.
A second method known as reputation was proposed by Wall and Mandel.92 In this
the last segment of a chain is removed and added to the other end. It has a faster
relaxation time than kink-jump methods but is not appropriate for dynamic studies.
A method that is good for very dilute systems is the pivot algorithm.87 It can be
used with very long chains to calculate the Flory exponent to high precision. In this
method a random site along the chain is tagged as the pivot point and one side of the
chain is rotated through a random angle.
The bond fluctuation method is also a well-used lattice algorithm for polymers.93
This uses a more complex lattice scheme with a coordination number larger than 100 in
order to capture some of the angular and stretching flexibility of polymers. It has been
used in a number of dynamic studies.94 Each move consists of a random segment
displacement (with a small maximum) and therefore resembles a Brownian dynamics
simulation.
Configurational Bias Monte Carlo
The method we use is a combination of Rosenbluth sampling with the Metropolis
method called configurational bias.95 It has a higher efficiency in simulations of long,
multiple chain systems than most other algorithms. This method will be discussed at
length in the next chapter as well as some extensions to it. The algorithm is as follows:
1. Randomly select a chain and a unit of that chain.
2. Erase all segments between that unit and the end of the chain.
3. Calculate the Rosenbluth weight (see below) Wold of the erased portion of chain.
4. Regrow the chain to its full length by randomly choosing directions that are
unoccupied.
5. Calculate the Rosenbluth weight of this new portion of chain Wnew.
6. Accept or reject this new portion of chain with probability Wnew Wold (or 1).
7. Go to step 1.
Where the Rosenbluth weight is designed to ensure the proper importance
sampling. It is calculated from the number of contacts with the surface and the solvent
that the section of chain has. For example, a segment touching the surface that has two
solvent contacts will contribute a factor: exp[-(Xs+2X)/kT].
This method is an improvement over those which use local moves such as
reputation, bond rotation, etc. Large sections of chain are moved at once, which decreases
level of "freezing" that occurs with the local moves. This allows the simulation of very
long chains.
We will discuss in Chapter 5, a significant advantage gained by a small
modification to the above algorithm. In step 1 we impose a maximum on the number of
segments that can be erased. This is useful for the simulation of long chains or dense
systems, since the probability of being able to regrow a chain that was cut in the middle
becomes negligible and is therefore a waste of time. In fact, the success ratio for
regrowths of a given size decays exponentially with the increasing size of the regrowth.
Also, instead of regrowing the chain from the original spot, it is regrown from the
opposite end from where the cut takes place. That way, the segments in the middle of
long chains will eventually have the opportunity to move unlike the standard method, in
which the interior segments of long chains are frozen.
Simulation of Bulk Properties
Bulk properties typically have faster relaxation times and are simpler for other
reasons as well. An excellent check of a simulation method is to test it against some of
the well known properties of bulk solutions. Critical exponents can be calculated such as
the Flory exponent (v) or the "susceptibility" exponent y. Here we will examine some
simulations that capture the correct scaling of the radius of gyration and end-to-end
distance.
Table II lists simulations performed on a single chain in good solvent. The size of
the simulation box must be several times the expected radius of gyration. Preferably 4-6
since fluctuations are of order of the rms radius of gyration. That poses a challenge on
the computer's memory for large chains. The Rg and Rend given for N=1-3 were
calculated by hand merely to show the limiting behavior at short chain lengths. Lattice
effects are obviously very important for short chains but appear only in the prefactor of
the scaling law at large chain lengths. The scaling laws Rg=0.41N.594 and Rend-0.99N0594
were found for the large chains (N=200-1000).
Table II. Properties calculated from the simulation of an isolated chain in an athermal
solvent.
N R g R end
1 0 0
2 0.5 1
3 0.6966 1.5492
5 0.9862 2.2760
10 1.5393 3.6579 0.6545
20 2.3739 5.7096 0.3569
50 4.1586 10.0402 0.1660
100 6.3050 15.1859 0.0952
200 9.5639 22.9356 0.0546
400 14.4466 34.6935 0.0317
1000 24.86975 59.59039 0.0155
The plot of this was shown in the Chapter 1 against corresponding simulations in
a near 0-solvent.. For the near 0-solvent the fitted power laws were Rg=0.57N.485 and
Rend=1.34N0.481. The power law develops even at chain lengths as low as 20 and 50, and
is accurate to less than a percent above 100.
Not only do simulations like this give confidence in the technique, but they
provide valuable information. The radius of gyration and end-to-end distance are
important length scales in all polymer problems and, as we will see, help to interpret
them. The overlap volume fraction (see Chapter 2) is also listed in the table and is a
critical parameter where many property transitions occur in polymer solutions. It is
N
calculated directly from the radius of gyration and chain length from 0* =
47R3 /3
Also, the prefactors of the power laws gives us the approximate correspondence
that the lattice size has to the Kuhn length of a polymer. For the 0-solvent an exact
correspondence of 1K = 1.48a is found by comparing the power law for the end-to-end
distance with the freely jointed chain. The end-to-end distance in a good solvent has a
prefactor of the same order of magnitude but slightly lower indicating that it describes
slightly smaller length scales.
To our knowledge simulations attempting to obtain the scaling of isolated chains
for N>1000 have not been done. While we did not use sophisticated extrapolation
procedures, our answer tends more toward the value 0.592 predicted from critical states
theory.11 At one time it was a great challenge to calculate the critical exponents with
simulations. Researchers extrapolated enumeration and Monte Carlo calculations to find
exponents for all dimensions and lattice types.21 Naturally it was found that only the
dimensionality of the lattice matters. The same critical exponents were obtained from all
3 dimensional lattices such as the simple cubic, diamond, and face centered cubic lattices.
Simulation of Adsorption
We have had success in the simulation of chains as long as 10,000 segments (see
the snapshot in Chapter 2). Most simulations reported for similar systems involve chains
only 50-100 segments in length. With simulation of multiple chains in equilibrium with
an adsorptive surface we are limited (for several reasons) to 1000-2000 segment chains.
Fortunately to obtain macroscopic properties of a polymer solution we do not
have to simulate 1023 (1 mole) chains. Simulations of bulk solution are done using a
small number (typically 10 100) of chains in a box with the usual periodic boundary
conditions that allow the chains to leave one side of the box and come in on the opposite
side. The box should be several times the average size of the chains (Rg). A box that is at
least 4 Rg wide will ensure, even in highly fluctuating systems, that the chains do not
interact with themselves through the periodic boundaries. As can be deduced from Table
II and the simulations described in Chapter 7, the smallest simulations boxes we used
were 3 4 Rg in the more concentrated systems (radius of gyration is smaller for these
systems) for N=200 and 1000 while the more typical case was 5 6 Rg.
We study the case of equilibrium between a bulk solution and adsorption onto a
flat surface. This is arranged by having periodic boundaries in the x and y directions and
2 impenetrable walls at z = 0 and L. L should be large enough so that the region in the
middle (z = L /2) is effectively a bulk solution. The symmetry of this situation is more
preferable to simulation with only one surface. There it would not be easy to create bulk
conditions even with a reflective boundary.
Equilibrium conditions. The basic criterion for equilibrium in bulk simulations is
that the simulation time be much longer than the relaxation time of the radius of gyration.
For simulation of adsorption, equilibrium conditions are harder to define or reach. The
characteristic relaxation times of surface conformations are typically longer than bulk
conformations.96 The main criterion we use however is that there be sufficient exchange
of chains between the bulk and the surface. This criterion becomes demanding for high
chain length in dilute solution since it literally requires 100's of kT to "rip" a chain from
the surface. At some separation distances there is more than one datum. The lower ones
are for simulations where the compressing surface is repulsive. A regular simulation can
never escape from a well of 100 kT. We have had difficulty in simulating chains of 200
or 300 segments using the traditional methods. When such chains have 50 or more
segments coupled to the surface, they are frozen to the local moves that traditional
methods are based on. It is the configurational bias methodology that can "tunnel"
through such steep activation barriers.
Solvent conditions. Athermal conditions can be ensured simply by setting the
monomer-monomer interaction energy, e, to zero. It is the simplest situation to study in
simulations and also it is well within the range of good solvent conditions (the most
relevant experimentally). Simulations are more naturally suited to study the effect of
excluded volume than theories. Recall the radius of gyration and end-to-end distances in
plotted in Chapter 2 which were obtained from simulations in athermal conditions. It was
much easier to find the Flory exponent than the exponent for a random walk.
That leads us to the question of 0-solvent conditions. A 0-solvent is defined as
the solvent that has a Flory-Huggins X= 1/2 which can be achieved experimentally by
varying the temperature of the solvent. For large or infinite length chains, this condition
results in Gaussian statistics. Thus in practice we can approximate 0-solvent condition by
varying the interaction energy between segments, e, so that the radius of gyration follows
the scaling law of an ideal chain, N1 2. Our best estimate comes from simulations for e =
0.25, 0.27 and 0.28 kT, yielding the following scaling laws at high chain length: Rend =
1.06N.53, 1.16NV.53, and 1.34N053 respectively. For these fits, the chain lengths 200, 400,
1000 were used for 0.25 and 0.27 and 400, 1000, 2000 for 0.28. A 2nd order polynomial
interpolation yields the 0-solvent condition: e=0.274 with Rend = 1.22NV500.
There could be subtle differences between a 0-solvent and a true random walk.
While e = 0.274 produces the square root law in the bulk solution, the prefactor may be a
function of concentration. As an example of the consequences, it is known that the
random walks confined in a slit or pore, will not expand laterally under compression.2
However, there is some evidence from simulations that finite length chains in a 0-solvent
do expand laterally.97
Markov matching: estimating the size of the lattice spacing a. In general an nth
order Markov process can be "matched" with a zeroth order Markov process. Alluded in
Chapter 2, this "Markov matching" (first proposed for polymers by Kuhn9) can give us
important information about the size of the lattice unit the chains on the lattice in
comparison to real polymer chains. Thus the quantities a / lc and M / N will be of
interest, where the former is the size ratio of the lattice to monomer unit and the latter is
the ratio of molecular weight to number of units in a lattice chain. In advance we know a
/ lc will be larger than unity because real polymers have more stiffness than lattice chains.
That is, polymer chains are a higher order Markov process.
This exercise consists of two parts: i) Matching the simulation chain to the freely
jointed chain using the correlations we just found for Rend in a 0-solvent; and ii) Matching
the chemical characteristics of real polymer to the freely jointed chain using either
experimental observations or statistical mechanics incorporating the chemical specifics of
various polymer types. In both cases, the matching proceeds as follows. Given N units
of an nth order Markov chain each of length 1, we match them to the 0th order chain (freely
jointed) which will have a smaller number NK units and a larger length per unit IK. The
conditions that the contour length L and the rms end-to-end distance Rend do not change
must be met.
L =NK=IN
(57)
= 1 /1/2 = clN1/2
Rend LKK 1 2
which leads to lK/I = c2 where c is the constant found either experimentally or in
simulation.
For case i) / is the lattice spacing a and we have found for isolated chains in a
0=solvent that c=1.22. This yields for the number of lattice segments in a Kuhn segment
lK/a = 1.5. This seems like a very reasonable answer since one would expect directional
memory of links to disappear somewhere between 1 and 2 links in a random walk on a
lattice. Bitsanis and ten Brinke98 have found for the same simulation method but in a
melt (which is also ideal) that the ratio is 1.1.
For case ii) we can use experimental measurements of Rg which have been
collated in the book by Fleer et al.4 Another option is to use calculations of what Flory
calls C_ = c2 which can be found in his second book.12 We will pursue the former now.
For a polymer with molecular weight M, where the molecular weight of each monomer is
m, we can use N = M / m in Equation 57. The length I = lc (chemical segment length),
will be obtained from the planar zig-zag conformation of polymers, using tetrahedral
valence angles of 109.50. The covalent bond length of carbon-carbon atoms is taken as
0.134 nm. Thus for a monomer that adds two carbons to the backbone of the polymer, lc
= 2(0.134) sin(109.5 / 2) = 0.2 nm.
I I
C C
C C
Figure 9: Length between carbons along the chain backbone.
Table III shows the results of these calculations for three types of polymer. In
general, the calculations by Flory for C are roughly the same as the experimental results
for k / lc. An improvement on the valence angle we chose could improve the agreement.
There are in these cases, between two and eight chemical segments in one our lattice
segments, and the level of detail of our lattice is approximately 1 nm. A 105 molecular
weight chain of polystyrene would correspond roughly to a 100 200 segment lattice
chain, while the same molecular weight version of polyethylene-oxide would be only 10
20 segment lattice chain.
Table III. Comparison of the simple cubic lattice to real polymer chains. Kuhn length, IK,
of the simple cubic lattice is IK = 1.49a.
R / Ml/2
Polymer/ Monomer c t R,/ Ik /c a a M/N M M
solvent weight (nm) (nm) N=100 N=200
PS / CH 104 10.2 0.0288 10.81 7.26 1.6 755.0 75501 151001
PEO / aq. 44 4 0.0343 2.88 1.94 0.6 85.2 8519 17039
PMMA/ (var.) 100 6.9 0.0261 8.53 5.73 1.3 573.3 57330 114660
PE/decanol 28 0.0435 6.64 4.46 1.0 124.9 12485 24970
PM/(var.) 14 6.7 3.32 2.23 0.5 31.2 3121 6243
PDMS / toluene 74 6 0.025 5.79 3.89 0.9 288.0 28803 57607
t C = 0 / NI2 calculated from rotational isomeric state model by Flory.8
tExperimental measurements collated by Fleer et al.4
Initial conditions. Before the attraction at the surfaces was turned on, the chains
were randomly grown. During the growth the chains also were subjected to local
relaxation moves (i.e., reputation, kink jump, rotation, and crankshaft) as outlined by
previously.98 These local relaxation moves have been shown to be effective for
simultaneous growth and equilibration, which is important for chains longer than 30 50
segments. The objective here is only to generate an initial configuration while the main
part of the system relaxation and collection of statistics are done with configurational
bias.
Polymer induced forces. There are two general methods99' 100 used to calculate
the equilibrium forces between a many body system and an object from a dynamic Monte
Carlo simulation on a lattice. We will discuss the latter, which is more accurate. Force is
related to the gradient of a potential energy field as f = -V Therefore, in a canonical
ensemble, the force required to maintain two flat surfaces at a separation H is related to
the derivative of the Helmholtz free energy, F, as
f(H) = F -F(H+a/2)-F(H-a/l2) -kTn Q+
dH a a Q_
We know that the canonical partition function, Q = exp(-U, /kT), can only be
obtained from simulation of systems for which the set of distinguishable configurations is
fully enumerable (note: there is a procedure that can estimate the partition function from
static Monte Carlo89). However, there is no fundamental reason why the ratio of the
partition functions of two similar systems cannot be obtained.
In practice it is not difficult to calculate the partition function ratio for two
systems of equal size, only differing in their potential fields. This means we cannot
actually move the wall from H + a / 2 to H a / 2. What is equivalent the approach
towards an infinite repulsive potential V in the layer next to the compressing wall. Thus
we wish to calculate the limit: Q(H; V = 0)/ Q(H; V- ).
Bennet101 developed a procedure for optimizing the accuracy in the calculation of
such a ratio (for an arbitrary system) through the identification of an adjustable weighting
function. Jimenez and Rajagopalan100 used this to obtain an implicit expression for the
change in free energy, AF, between two lattice systems with difference in potential at the
upper wall of AV
kmax pb a (AF-kAV)/kT
k (AF-kA )kT =0. (59)
k=0 + e- /
All that is required to solve Equation 59 for the change in free energy between
two systems is the probability distribution of surface contacts for each system: Pk. This is
just the probability that there are k segments at the compressing wall. Equation 59 is
valid for any change in potential energy AV If this change is infinite, the equations
reduce to
AF/kT = InPa
(60)
and all that is required is one simulation for V = 0 to find the probability of zero contacts.
In practice, if the probability of zero contacts is very low (i.e., strong force dense
systems) then it cannot be obtained accurately. In that case one incrementally changes
the repulsive potential V and calculate AF at each step from Equation 59. Eventually, for
some highly repulsive potential, the probability of zero contacts can be found and thus
the final term from Equation 60 is added.
We summarize calculation of force at a distance H as follows:
1. Simulate the system with a certain potential at the upper wall.
2. Simulate the system again with a slightly higher potential. The change in potential
must be small so that the new distribution overlaps with the previous one allowing the
solution of Equation 59 for the incremental change in free energy.
3. Repeat 2 until the repulsive potential is high enough that the probability of zero
contacts can be calculated with accuracy. The final increment in free energy comes
from Equation 60.
The simulation of two physisorbed layers is a little more difficult since one must
start from a potential at the upper wall that is attractive, and increment gradually in the
repulsive direction until extrapolation to infinity is possible. This is only stage one
however since we have only calculated the free energy difference between a system with
an adsorptive upper wall at z = H + a / 2 and a system with a non-adsorptive upper wall
located at z = H a / 2. In stage 2 we need to find the free energy difference of re-
adsorbing the chains. This is done by incrementally decreasing the potential, starting
from zero, and continuing until the original adsorption strength is obtained. For dual
physisorbed layers, the energies from stage 1 and 2 are large number with opposite signs
81
which is the reason for poor accuracy for the calculation of weak forces ( < 0.001 kT/ a3)
when the separation is large. The number of simulations this requires altogether can be
as many as 15 or 20. Nevertheless, we have had some success in formulating this
technique to obtain unique information on bridging forces.102
CHAPTER 5
REPTATING CONFIGURATIONAL BIAS MONTE CARLO
FOR LONG CHAIN MOLECULES AND DENSE SYSTEMS
Introduction
The exponential attrition of configurational bias for long chains is reduced to
quadratic by a simple modification of the basic move. Each trial move begins instead on
the side of the remaining sub chain opposite to the cut location. This type of move is
akin to a large-scale reputation and in the limit of a one-segment cut is reputation. The
extension is proven to require the same Rosenbluth weighting scheme as the original
algorithm. Several examples are used to demonstrate the improved performance of the
new method. A single chain is analyzed in isolated and more concentrated environments.
Also, some characteristics of the more demanding problem of polymer physisorption are
elucidated. Finally, a further extension to the method that forces repeated selection from
one end of the chain only is considered. The method while having the advantage of
completely eliminating attrition is rigorously incorrect. However it may be teleologically
correct in many situations.
Development of algorithms for sampling statistics of long chain molecules has
been a challenge since the first days of Monte Carlo. One of the most successful models
because of its simplicity is a self-avoiding walk on a cubic lattice. While this model
suffers in its description of microscopic properties it has been used to analyze large-scale
aspects of many important problems. It corresponds to a coarse-grained polymer chain in
an athermal solvent. Sometimes nearest-neighbor contact interactions are used as well.
This problem will be the subject of much of this chapter, although most of our results can
be applied to off-lattice models as well.
Many of the early methods for polymer chains used static Monte Carlo variants in
which each microstate of a system is generated independently. The first method used for
polymer was called "simple sampling." Because Boltzmann statistics require non-
overlapping and overlapping configurations to be generated with equal probability the
method fails for all but exceptionally short chains. The success rate for chain creation
faces a dramatic exponential attrition with increasing chain length.
"Enrichment" and "biased sampling" were among the more successful attempts to
overcome this problem. However their effect is only to "postpone" the inevitable
attrition. Biased sampling, developed by Rosenbluth and Rosenbluth,89 is the basis of
configurational bias. In contrast to simple sampling, biased sampling allows a choice
only from unoccupied neighboring sites for the sequential placement of trial chain
segments. These early lattice methods are reviewed at length by Kremer and Binder.87
Ultimately complex polymer systems yield their secrets only to dynamic (i.e.,
Metropolis) Monte Carlo.103 "Pivot," "kink jump" and reputation algorithms have all had
some degree of success in simulating longer chains and denser systems and can even be
combined. As opposed to these "traditional" Metropolis methods, configurational bias is
a more "forward-looking" algorithm. It determines each trial configuration with respect
to the system's current configuration in a way that hopefully produces a microstate with
frequency more proportional to its weight in the Boltzmann distribution of microstates.
Developed originally for a canonical ensemble on a lattice104' 95 it has been extended to
off-lattice,105, 106 grand canonical ensemble and Gibbs ensemble systems.107 These
techniques have been reviewed in the book by Frenkel and Smit,5 as well as a
generalization the configurational bias method to all (i.e., not necessarily polymeric)
molecular systems.
Several advances have appeared since then in CBMC. Algorithms which form
trial moves by regrowth of internal instead of end segments have received much
attention.5, 108, 109, 110, 111 The first implementation was on a lattice by Dijkstra et al.108 and
uses a priori knowledge of the closure probability from the properties of a random walk.
This class of algorithms have the potential to generalize CBMC to polymer systems of
many topologies. Applications exist towards simulation of long linear chains, branched
chains, long grafted chains and ring chains to name several problems that are presently
difficult because of their low concentration of endpoints. However the growth of a sub
chain between two fixed points forces a sacrifice of the method's simplicity. Another
advanced extension that overcomes some of CBMC's inherent limitations is the recoil
growth method.112 This is useful in dense, long-chain systems since each growth in the
trial move uses a "retractable feeler" to look several steps forward as compared to the
one-step pre-cognizance of conventional CBMC.
CBMC has been implemented in many phase equilibrium studies.113 Therefore
serious attempts to improve its efficiency for large systems with complicated interaction
potentials have been made.114' 115, 116, 117 Interactions between long chains and interfaces
have been largely ignored by CBMC and require specialized methods.
The organization of this chapter is as follows. In the next section two extensions
to CBMC are proposed and analyzed theoretically. These make long chain simulations
easier through very simple means. The following section contains several comparisons
between conventional and new CBMC in terms of accuracy of results and performance.
This is succeeded by a section summarizing our main findings.
Proposal and Detailed Balance of Configurational Bias Extensions
End-Switching Configurational Bias
In this communication, certain issues are addressed regarding modifications of the
original configuration bias Monte Carlo (CBMC) algorithm proposed by Siepmann and
Frenkel.95 We analyze two modifications that are aimed at relaxing more effectively the
interior segments of long chains in dense systems. In general these methods involve the
removal and random regrowth of pieces of the chains in the system. In contrast to the
original CBMC algorithm (Figure 10a) we propose to grow the trial sub chain opposite
the end from which a sub chain is cut, as shown in Figure 10b. The advantage of this
method is in the relaxation of the interior segments of the chain. Reptation dynamics tells
us that the rate of backbone renewal will increase quadratically with the number of
segments in the chain instead of exponentially. That is to say, the average number of
successful end-switching CB moves that produce a complete renewal of the chains
backbone is proportional to the square of the number of subchains, -(N/Nc)2 N2,
which amounts to a much milder attrition.
For all Monte Carlo methods that generate ensemble snapshots sequentially, the
primary criterion is that the average transition rate between any two configurations must
be equal in the forward and reverse directions. This condition is known as "detailed
balance." A microstate (or configuration), I say, is defined exclusively by a set of
n(N +1) position vectors, where n is the number of chains and N is the number of links
in each chain. For simplicity, lets denote any such set by the notation rf The number
density, PI, of such configurations in the equilibrium ensemble is proportional to the
Boltzmann factor, exp[-p&D(r,")]. As usual, P and ( refer to the inverse thermal energy
and the potential energy respectively. In what follows, we will consider only athermal
systems in which the Boltzmann factor is either 0 if overlaps are present in the system, or
1.0 if no overlaps exist. We can consider the transition rates between arbitrary states I and
J as the product of the number density of each state and the probability, H, of reaching J
from I, or vice versa, via the simulation algorithm. Therefore detailed balance is written
as
PT1,, = PHn,,. (61)
This transition probability can be separated into two factors due to the probability, TI1, of
generating trial state J from state I and the acceptance / rejection rule, AI1 for this trial
move. Normally Monte Carlo methods are devised by proposing a set of moves and thus
fixing TIZ, and then finding an acceptance rule that is consistent with Equation 54.
It is useful to consider the example of two 4-mer configurations shown in Figure
11. For simplicity we use two-dimensional chains inscribed on a square lattice.
Additionally, we will only incorporate removal / regrowth moves of Nc = 2 links at a
time (normally Nc is picked at random). What is the probability, TI1, for for generating
trial move J from I? First a chain in the system is selected at random, and then one end of
the marked chain is selected at random. Thus a factor (1/2)(1/2) is contributed towards
TIj in the example since there are two chains (and of course, two chain ends). At each
stage of the trial insertion, J, a link is added randomly to one of the z available directions.
Thus TIj = (1/2)(1/2)(1/3)(1/3) for the forward move and for the reverse move, Tj, =
(1/2)(1/2)(1/2)(1/2). Therefore the acceptance rule must compensate for the fact that
moves from J to I in this case are more likely. Specifically, A1I / Ajj, is required to be
[zz 2]J /[Z z2] = (2- 2)/(3- 3) 4 / 9.
A generalization for the acceptance rule to athermal lattice systems of any size is
suggested by the above example. An arbitrary configuration J is generated from I with a
probability,
Tj const (62)
n2[Z1 Z2ZN1 W,
where Wj is the Rosenbluth weight generated by the trial portion of the chain as it is
grown in configuration J. That is to say, Wj = W(Nc ) where W(i +1) = (zI / 5)W(i) and
the factor of 5 corresponds to the number of non-backfolding directions on a cubic lattice
(for a square 2D lattice it would be 3). Therefore we find the acceptance rule, in the
absence of any energetic interactions, must satisfy,
Tjl WJ A,
(63)
TI, W, A,
and the best choice for the acceptance rule is therefore
A,j = min[1, W, / W,] (64)
This is exactly the expression used in standard CBMC (without end-switching) and thus
constitutes a proof that they are the same in terms of weighting functions.
We should point out that in these algorithms, the weighting functions a priori
have Nc factors in them, not N factors. That the N Nc factors representing the uncut
portion of the chain, if included into the weighting functions, would cancel out in the
standard method but would not this end-switching method, is purely coincidental.
Persistent End-Switching Configurational Bias
One might entertain the notion that by continually cutting segments from the same
side of each chain in a system. An advantage would be gained in the relaxation of
exceptionally large chains since a chain renewal would be achieved in exactly (N/Nc)
moves. It is tempting to raise the objection that, due to the fact that immediate reversal of
a move is impossible, the scheme automatically violates the law of detailed balance.
However in reality the principle of detailed balance makes no such implication. This
misconception might be due to the popular usage of the term "microscopic reversibility"
when the term "detailed balance" is more a realistic description. One example in previous
literature of a scheme that is "non-reversible" is the persistent reputation algorithm118
(however they switch directions when the chain meets a dead end). As stated above,
detailed balance is only a sufficient condition. The necessary and sufficient condition is
sometimes referred to as "semi-detailed balance" which is less restrictive but harder to
prove.105 This just states that the exchange between two states I and J must be equal.
The difference between detailed balance and semi-detailed balance can be non trivial
(and several examples exist to demonstrate this).
It is best to clarify the persistent end-switching CBMC proposal by the use of an
example. Figure 12 shows two possible configurations of a system of tetramers, again on
a two-dimensional lattice. The trial move algorithm is to remove links 1 and 2 and add
them randomly to link 4. If such a move is accepted, the asterisk (denoting head of chain)
is updated. Any state I in the ensemble then is either of two kinds as regards asterisk
position, which we will call sub states i and i'. For property calculations we do not need
to distinguish between these two sub states, but there certainly is a distinction in terms of
the transition probabilities generated by the algorithm of choice. In Figure 12 the
probability, T,, is (1/2)(1/3)(1/3) and the probability T7,,, is (1/2)(1/2)(1/2). All of the
other combinations, i.e.,, T7,,, Tj,,, Tj,,, T, are zero since these transformations cannot
occur. For instance, i cannot be reached from j' by our proposed move. The only
possibilities are a transition from i toj or a transition from' to i'.
For an arbitrary athermal lattice system it is necessary to recognize that the
requirement of detailed balance still holds in the form of Equation 54. A detailed balance
on sub states is not required! If it were, then the departure of this algorithm from the
detailed balance condition would be of the most extreme kind. All that is required is that
the "flow rate" between states is balanced in the reverse and forward directions. We have
shown that the flow rate in the forward direction consists only of i toj transitions. Similar
is the reverse direction flow rate, which completes the detailed balance expression,
P'ni, = PnH (65)
The obvious question is whether the number of sub states are equal in any given
state, i.e., does P, = P,= PI? If this were true then in the absence of energetic
interactions all configurations are equally likely. Consequently P, = P and an expression
such as Equation 63 would define the appropriate acceptance rule, and the same weights
as in Equation 55 would produce a correct Monte Carlo algorithm. However this
assumption is unfounded. The condition can only be enforced either by a detailed or
semi-detailed balance on sub states; clearly it violates the former condition. In what
follows we show that the persistent end-switching scheme must even be in violation of
the semi-detailed balance condition. This is accomplished by picking illustrative
examples for which P, # P, rather than by a general proof.
Let us attempt to find an example microstate for which the algorithm a priori will
not generate equal numbers of chains with the asterisk located on either end. Again we
shall use only a two-dimensional square lattice to simplify matters. A first attempt is to
consider a fully extended dimer normally adjacent to a surface (or a line). This is shown
in Figure 13a with the asterisk located at the far end and in Figure 13b with the asterisk at
the near end of the chain. Clearly there are 3 possible configurations from which the
former can be generated and only 2 for the latter, and therefore the second configuration
should be 2/3 as likely as the first.
The second example, depicted in Figure 14, also shows that sub states are non
ergodic. The same dimer is now up against an object (either an irregular surface or
another chain) that traps its last segment. Therefore the configuration j with asterisk in
the trap cannot be reached by a single segment move. The opposite configuration j' can
be reached from 3 other possible configurations. If the trap is a fixed object than it is an
example of a non ergodic configuration and there will be zero such states. Even if the
trap is not fixed, i.e., it is another chain, the microstates I and J do not balance.
Tests of the Proposed Algorithms
The methods will first be assessed by the simulation of an isolated chain. The
accuracy of the end-to-end distance of the chain and the associated scaling exponent will
be examined, and in a more sensitive test, the segmental pair correlations and radial
density gradient will be analyzed. This is followed up with the analysis of two further
systems. To illustrate the effect of system density on the algorithms, a semidilute and
concentrated bulk solution will be analyzed. Finally, a problem which inspired much of
this work will be outlined. This is the physisorption of polymers onto a flat surface in
equilibrium with a bulk reservoir of polymer solution. To facilitate discussion the
abbreviations: CBO, CB1 and CB2 will often be used in reference to the standard, random
end-switching and persistent end-switching CBMC methods.
|
