journal of
MO LECU LA R
LIQUIDS ELSEVIER
Polyelectrolyte
Journal of Molecular Liquids, 61 (1994) 209-220
chain
dimensions
K. Krishnaswamy and M. Fixman Department of Chemistry, Colorado State University Fort Collins, CO 80523, USA
Abstract The parameters that enter excluded volume theory, namely the persistence length and the excluded volume between segments, depend in a qualitatively obvious way on the polyelectrolyte charge and the solution ionic strength. However, a quantitative understanding is quite difficult to achieve for realistic models of flexible chains. Therefore a simulation of these parameters was performed for small sections of polyelectrolyte chains and the results were incorporated into excluded volume theory. The theory was suitably modified to take into account the fact that a segment is a finite rather than an infinitesimal fraction of the whole chain, as is usually assumed. Excellent agreement was found between the predicted expansions of whole chain dimensions and simulation results.
i. I N T R O D U C T I O N For one very simple model of a polyeleetrolyte chain the dependence of the persistence length and excluded volume parameters on the backbone charge and ionic strength may be calculated by analytical or numerical techniques. This simple model, a gently curved and uniformly charged cylinder, has a rather limited application to real systems and it seems unlikely that purely theoretical techniques will succeed in the treatment of realistic models. The purpose of this work is to illustrate an alternative method of dealing with more complicated polyelectrolyte models by the computer simulation of small sections or segments of the chain. Such simulations can provide values for the parameters in excluded volume theory, namely the persistence length and the excluded volume between segments, and are feasible because these parameters reflect local rather than global properties of the chain. Therefore only small sections of the chain require simulation. Simulations of the whole chain have been performed in order to test the theory. Of course many simulations of polyelectro-
SSD10167-7322 (94) 00761-K
210
lyre chains have been carried out previously [i], but generally with the objective of testing theoretical treatments that were based on different and simpler models than were used for the simulations. The most useful current theory of polyelectrolyte chain expansion was presented many years ago by Odijk and Houwart [2] and is purely analytic. Odijk and Houwart based their theory on a worm model and considered separately the two sources of chain expansion in a "two parameter" theory [3]. These are (a), the increase in persistence length due to electrostatic repulsion between charges which are close together along the backbone, and (b), the excluded volume interactions between distant charges. The model itself and the restrictions imposed in order to allow an analytic treatment limit the application of the theory. In order to understand these limitations we briefly summarize the two ingredients of the Odijk-Houwart theory. First, the dependence of persistence length on backbone charge and ionic strength has been studied in detail for the uniformly charged worm model in both the Debye-HOckel [4,5] and the Poisson-Boltzmann approximations [6,7]. The Debye-HOckel approximation may be solved analytically subject to the limitation that p>>~-l>>a, where p is the persistence length, - I
is the Debye screening length, and a is the radius of the
wormlike chain. For uniformly charged curved cylinders or toroids the numerical solutions of the Poisson-Boltzmann equation indicate that "counterion condensation" [8,9] can be accounted for, to a reasonable approximation, through the use of an effective backbone charge. Second, the further expansion of the chain due to repulsion between charges far removed from each other along the backbone was treated by OdlJk and Houwart according to the ordinary theory of excluded volume effects for nonelectrolyte polymers [3,10]. This theory is basically a perturbation theory, possibly a renormallzed perturbation theory, in which any particular pair interaction between chain segments has a small effect on the overall expansion. In the two-parameter version of excluded volume theory used by Odijk and Houwart the chain must be sufficiently long that Gaussian statistics apply in the absence of excluded volume interactions. The combined theory has had broad success for polymers, such as very long or very short DNA, that resemble the model and satisfy the rather
211
s e v e r e c o n d i t i o n s imposed on c h a i n l e n g t h , p e r s i s t e n c e
length,
and Debye
s c r e e n i n g l e n g t h . A d e t a i l e d u n d e r s t a n d i n g o f o t h e r models o r i n t e r m e d i a t e chain lengths is still
quite limited,
and i s u n l i k e l y t o advance much
further on the basis of purely analytic treatments. However, the conceptual distinctions that underlie the Odijk-Houwart work seem reasonable for any model and will be used here. In outline our calculation proceeds as follows. First, divided up into sections labeled I, ISISN,
the chain is
each of which may contain many
monomer units at positions labeled, collectively,
{RI}. The total electro-
static free energy of the chain is assumed in this work to be a palrwlse sum of known energy terms WIj that depend only on the coordinates of sections I and J. The backbone model may be any rotational-isomeric state model. Electrostatic interactions are approximated as screened Coulomb interactions with effective charges. In principle a simulation of the ion atmosphere could be performed in order to improve on the screened Coulomb interaction, but this is not attempted here. Second, a reference model for the application of excluded volume perturbation theory is defined by a restriction of the free energy to self interactions WII and nearest neighbor interactions WI,I± I. If the chain has a sufficiently large number of sections then the distribution of end-end distances in the reference model becomes Gausslan. But regardless of the number of sections or the complexity of the full model, the restriction to local interactions in the reference model makes feasible a simulation of its conflgurational
statistics. Third, a first order excluded volume perturbation theory of those interactions discarded in the reference model is carried out. For chains sufficiently short that first order excluded volume theory suffices, a direct comparison between the theory and full chain simulations is performed. For longer chains we have to infer the approach of excluded volume parameters to their asymptotic values, in order to apply renormalized excluded volume theory. The various models and interaction energies that are necessary to define the reference chain and the full chain are introduced in Section 2. The modified perturbation theory is presented in Section 3, and the results and the particular choices of parameters for which they were obtained are
212
presented and discussed in Section 4. Our conclusion is that the combination of simulations of short sections of the chain with excluded volume perturbation theory is capable of providing a theoretical understanding of conformational statistics for a wide variety of polyelectrolyte models. There still remains the need to incorporate an explicit simulation of the ion atmosphere in order to avoid the use of effective backbone charges.
2. C h a i n M o d e l
and Computational
Algorithms
The backbone model consists of a rotational-lsomerlc chain of n+l beads. In practice the bond lengths and bond angles were constrained to uniform values b
and 0, respectively, and rotational potentials for the o dihedral angles were omitted, but these restrictions are not essential to the approach. Each bead i is assigned a charge qilel, where e is the electron charge. The interaction between beads was restricted to a screened Coulomb potential uij (Rij). Dimensionless variables are used in which kBT is the unit of energy and b
o
is the unit of length. Then
(2.1)
uij (Rij) - Bqiq j exp(-mRij )/rij,
where b B is the Bjerrum length, i.e., o B - e2/(kBTDwbo ) .
(2.2)
At T-298K and with the dielectric constant of water chosen to be D -78.5, w
b B-7.136A. The q's have been treated as effective charges which should be o
chosen somewhat smaller than realistic values in order to account for counterion condensation. However, a considerable spread of values has been examined in order to examine the behavior of the theoretical approach. In this work b -IA throughout. Although it would be possible to scale o
the charges and Debye length to reflect another value, our focus has been on tests of the approach rather than on comparisons with experiment. More fundamental than the particular choice for the backbone model is
213
the distinction between the reference chain model and the full chain model. This distinction underlies the perturbation theory. The definition of the reference model is completed as follows. Each consecutive sequence of K bonds constitutes a section. The last bead of a section is a "phantom" bead because it is uncharged and because its physical position in the full chain coincides with the first bead of the next section. The first bead of a section, rather than the last (phantom) bead of the preceding section, carries any physical charge that may be present. The sections are numbered I, i ~ I ~ , where Nmn/K is required to be an integer. In the reference model the interactions are restricted to self-interactions beween the beads of each section and neighbor interactions beween the beads of adjacent sections. The full chain model includes all interactions. Compared to the full chain, the reference model has relatively simple statistics due to the restriction to nearest neighbor interactions between sections. Sufficiently long reference chains will have a G a u s s i a n distribution of the end-end vector and may be expected to behave qualitatively llke the worm model, especially for strong Coulomb repulsions which favor a parallel arrangement of adjacent sections. Several restrictions on the choice of K must be observed. The root mean square length I / 2 ^ of section I, where hi is the vector from the first bead to the last (phantom) bead of section I, must be comparable to or somewhat greater than the Debye length in order to minimize interactions between sections I and I±2. Because all interactions between segments I and J for ll-Jl~2 will be treated by the perturbative
methods of polymer excluded volume theory, it
is necessary that the total pair interaction between I and J be small or, for large ll-Jl, that it occurs with low probability. Except for that limitation smaller sections are preferred because their conformational distribution must be inferred from simulations rather than analytically. A section may be much smaller than the persistence length; the latter can grow without bound as the backbone charges are increased to arbitrarily large values. Simulations of either the reference chain or the full chain were based on the Metropolis algorithm. Configurations of the reference chain were generated by two alternative methods, either by simulation of a long
214 chain with interactions restricted to nearest neighbor sections, or by simulation of individual sections and subsequent coupling of these sections, with proper coupling energies and weights, to form instances of a long reference chain. The latter method requires umbrella sampling if adjacent sections have strong electrostatic interactions. Initial configurations were chosen to be either those for a random distribution of dihedral angles or, for checking, fully elongated chains. In either case the configurational phase space of the chain is restricted to its dihedral angles and standard methods [ii] allow the spatial configuration to be determined from these angles.
3. M O D I F I E D
PERTURBATION
THEORY
In the usual perturbation theory of polymer excluded volume effects [3] the chains are assumed to be very long and the conformational statistics of the reference model are assumed to be Gaussian. Since the theory will be tested here by comparison with simulations of rather short chains, these simplifications cannot be used. In other respects the perturbation theory is standard and a brief exposition is considered sufficient. Let WIj be the screened Coulomb interaction energy of all beads in section I with those of section J, i.e. WIj is the sum of uij over beads i belonging to section I and beads j belonging to section J. Since terms WII and WI,I+ 1 are included in the reference model, the probability density P{R} of the full chain is related to the probability density Po{R} of the reference model by
P{R} = Po{R}exp(-AU)/o,
(3.1)
where <>o indicates an average taken with weight Po{R} and AU is given by
AU - ~'~Wlj.
(3.2)
The notation ~'~ here and below indicates a double sum of I and J over the N sections, subject to the restriction (J-l)~2. With the cluster function flJ defined by
215
flJ " exp(-WIj)-I'
(3.3)
Eq. (3.1) may be expanded through "single contact" terms as
•
f
(3.4)
.
We have found it useful for a few extreme choices of parameters to make o-0 by adding a constant AIj to the WIj functions for small II-Jl. AIj depends on I and J but not on the configuration,
and it is apparent from
Eqs. (3.1) and (3.2) that the equations are unaffected by this change of reference energy. However, the change improves the convergence of the perturbation expansion when WIj contains large terms which fluctuate only slightly. This situation occurs for large charges, large Debye length, and short sections. Equation (3.4) has been applied to the calculation of the mean square end-end distance
-
o + Z'EIo-oo ] + ...
(3.5)
and also to the perturbation expansion of the characteristic function of the end-end vector. We have concentrated on in this report because excluded volume perturbation theory is most extensively developed for that quantity. For any chain length let the excluded volume expansion parameter
2 n
be defined by
2 n
R2
" < >I<
R2
(3 6)
>o '
In the long chain limit the ratio 2 ~ 2 ( z ) , 2 n (z) of the reduced parameter z defined by z = [3/(2~(R2>o)]3/2~N 2,
i.e.,
2 becomes a function n
(3.7)
216 where ~ is the effective volume excluded from one section by the presence of another section and has the form,
- ~(l-exp(-W)>odR.
(3.8)
Here W is any of the WIj functions all of which are assumed to have the same form. The integral ~ is one of the factors in -o evaluated for very large
a
ll-J[. The expansion of a 2 starts with the terms 2
[3]
- i + (4/35z + ....
and Muthukumar and Nickel
(3.95
[i0] have constructed an expression for a 2 based
on a series analysis of their sixth order perturbation calculation,
ln(= 2) - O.17722n(l+7.524z+ll.O6z25.
(3.10)
Equation (3.10) and indeed the dependence of a on the single parameter z are restricted to chains of sufficient length that not only ° but every significant
interaction in the reference model can be evaluated on the
basis of Gaussian statistics. z to N I/2.
In that limit
We have adopted a semi-emplrical
is proportional
to N and
o scheme in order to apply excluded
volume theory to chains of modest length. We start with the actual first order perturbation result for finite chains, which is known from Eq. (3.5) 2 and the simulation of segments. The expansion factor u may be expressed n as
2 an
m
I + (4/3)z I + ....
which constitutes,
along with Eqs.
(3.11)
(3.5) and (3.6), a definition of z I. The
asymptotic form of z is known for long chains, namely z-cN I/2, and the parameter
c in that equation is also known from Eqs.
the simulation of uncoupled segments.
(3.75 and (3.85 and
It turns out of course that Zl
mainly because Gaussian statistics for the contacts between segments I and
217
J is not realized until
[l-J I is very large. For small values of [l-J I the
probability of contact is much less than given by the Gaussian estimate. An extensive analytic study of the relationship between z I and z for the worm model was performed long ago by Yamakawa and Stockmayer [12], but we have contented ourselves with the approximation
z I = cN I/2 - d,
(3.12)
with d independent of N. The constant term d in Eq. (3.12) is expected to be only a leading order correction, valid for large N, but in practice Eq. (3.12) is found to be a generally adequate approximation. It follows that substitution of z I for z in Eq. (3.10) will force long chain excluded volume theory to agree with the simulations of short chains. On the other hand the distinction between z I and z becomes unimportant for long chains where the expansion is dominated by interactions between segments I and J for which II-J I is large. In general, then, we expect that 2 2 n ~ (Zl)'
where the function ~
2
(3.13) .
is obtained by substitution of z I for z in Eq.
(3.10), will be an adequate interpolation formula for all chain lengths greater than the segment length. Tests of Eqs. (3.12) and (3.13) will be summarized in the next Section.
4. R e s u l t s For all the results to be reported the backbone model consisted of a chain for which bo-IA, 0-I09.5 ° , and an equal charge q was placed on every other bead, starting with the first. The sections of the reference model always consisted of I(-20 bonds/sectlon. Consider first the relationship between the excluded volume parameter z for a very long ("infinite") chain and the parameter z I for a finite chain. For this comparison z I was evaluated from the first order perturbation terms in Eq. (3.5) for a chain with q-0.2 in units of lel, ~-I-IA and
218
for n ranging from 120 to 480. The value of z was obtained from Eq.
(3.7),
with o taken to be the actual mean square end-end distance for the reference chain; ~ was obtained from a separate simulation of the interaction between two properly charged sections.
The results are shown in
Figure I, with filled circles for z and open circles for z I. As anticipated above, z-z I is approximately
independent of N. z I vanishes
in the vicinity
of N=3 for a chain with the specified structural and charge parameters.
Figure I. The asymptotic excluded volume parameter z and the first order parameter z I vs. JN, the square root of the number of segments.
1.2
,
,
,
,
,
OQ
1.0 Z
0.8 0.6
© O
0
0.4
0
Z1
0.2
0
0
©
00.
'
2.0
2.5
'
'
'
3.0
3.5
4.0
I
4.5
5.0
219
-gaff 2(Zl) with
Next we compare in Figure 2 the predicted values
the actual values of a 2 determined from a simulation of the whole chain n with all interactions included. The open circles represent the same data shown in Figure I. The triangles were all evaluated for a chain with n-120 bonds, ~-I-6A, and with q varying from 0.3 to 6 units of lel. The dotted line is just the result of first order perturbation theory, Eq. (3.11), and the solid line is the Muthukumar-Nickel expression, Eq. (3.10). Again, the agreement seems excellent.
Figure 2.
2 vs. the first order excluded volume parameter z I. n
2.0
I
I
I
I
l
l
1.8
0 0
1.6
1.4 1.2 1.0
~, 0.0
O. 1
I
I
I
I
I
0.2
0.3
0.4
0.5
0.6 Z1
0.7
220
5. C o n c l u d i n g
Remarks
This paper is dedicated to the memory of our friend and colleague, Teresa Fonseca.
6. R e f e r e n c e s 1. C. E. Reed and W. F. Reed, J. Chem. Phys. 97, 7766 (1992), and references therein. 2. T.Odijk and A.C.Houwaart, J.Polym. Sci. Polym. Phys.Ed., 16, 627 (1978) 3. H.Yamakawa, "Modern Theory of Polymer Solutions" (Harper and Row, New York, 1971). 4. T. Odijk, J. Polym. Sci. Polym. Phys. Ed. 15, 477 (1977). 5. J. Skolnlck and M. Fixman, Macromolecules 10, 944 (1977). M.Fixman and J.Skolnlck, Macromolecules 11, 863, 1978. 6. M. Le Bret, J. Chem. Phys. 76, 6243 (1982). 7. M. Fixman, J. Chem. Phys. 76, 6346 (1982). 8. G. S. Manning, J. Chem. Phys. 51, 924 (1969). 9. C. F. Anderson and M. T. Record, Jr., Ann. Rev. Phys. Chem. 33, 191 (1982). 10. M. Muthukumar and B.G.Nickel, J.Chem. Phys. 86, 460 (1987). ii. P. J. Flory, "Statistical Mechanics of Chain Molecules" (Hanser Publishers, New York, 1989). 12. H. Yamakawa and W. H. Stockmayer, J. Chem. Phys. 57, 2843 (1972).