J. Phys. Chem. B 2010, 114, 3185–3196
3185
Electrostatic Persistence Length Marshall Fixman Department of Chemistry, Colorado State UniVersity, Fort Collins, Colorado 80523 ReceiVed: January 20, 2010
The persistence length is calculated for polyelectrolyte chains with fixed bond lengths and bond angles (π θ), and a potential energy consisting of the screened Coulomb interaction between beads, potential wells Rφi2 for the dihedral angles φi, and coupling terms βφiφi(1. This model defines a librating chain that reduces in appropriate limits to the freely rotating or wormlike chains, it can accommodate local crumpling or extreme stiffness, and it is easy to simulate. A planar-quadratic (pq), analytic approximation is based on an expansion of the electrostatic energy in eigenfunctions of the quadratic form that describes the backbone energy, and on the assumption that the quadratic form not only is positive but also adequately confines the chain in an infinite phase space of dihedral angles to the physically unique part with all |φi| < π. The pq approximation is available under these weak constraints, but the simulations confirm its quantitative accuracy only under the expected condition that R is large, that is, for very stiff chains. Stiff chains can also be simulated with small R and small θ and compared to an OSF approximation suitably generalized to chains with finite rather than vanishing θ, and increasing agreement with OSF is found the smaller is θ. The two approximations, one becoming exact as R f ∞ with fixed θ, the other as θ f 0 with fixed R, are quantitatively similar in behavior, both giving a persistence length P ) P0 + aD2 for stiff chains, where D is the Debye length. However, the coefficient apq is about twice the value of aOSF. Under other conditions the simulations show that P may or not be linear in D2 at small or moderate D, depending on the magnitudes of R, β, θ, and the charge density but always becomes linear at large D. Even at a moderately low charge density, corresponding to fewer than 20% of the beads being charged, and with strong crumpling induced by large β, increasing D dissolves blobs and recovers a linear dependence of P on D2, although a lower power of D gives an adequate fit at moderate D. For the class of models considered, it is concluded that the only universal feature is the asymptotic linearity of P in D2, regardless of flexibility or stiffness. I. Introduction 1
2
Quite some time ago Odijk and Skolnick and Fixman (OSF) calculated the electrostatic contribution to the persistence length P of a charged, wormlike chain in salt solution. With the aid of several restrictive assumptions,3 OSF concluded that the increase in P due to Debye-Hu¨ckel (DH) repulsion between the charges is proportional to D2, where D is the Debye screening length. The OSF result was generally accepted4-6 as correct for stiff, wormlike chains, but it was of course understood that applications of the result to experimental measurements might reveal inadequacies in the model, such as failure of the DH approximation7-10 for highly charged chains, or failure of the structural assumptions.3,10,11 Following extended debate over the course of decades, a consensus5,12-14 seemed to emerge on the basis of simulations and variational theory that the OSF result was also correct for discrete models of flexible chains, to the extent that a quadratic power law dependence of P on D would occur at large D. This consensus has recently been expanded by work15,16 that supersedes earlier objections.17-20 The nature of the consensus and the evidence for it are by no means entirely satisfactory. Asymptotic power laws are only part of the story, and a consensus is not a proof. The theoretical models are too simple to expect quantitative agreement with experiment, or to regard agreement as other than accidental if found, and yet few theoretical models are designed to allow quantitative comparison with simulations. The present work is intended to address these deficiencies, especially the last.
In brief summary, potential energy functions are introduced into the dihedral angle space of freely rotating chains, first to cover the range of uncharged models from highly flexible and locally curled at one extreme, to highly elongated and stiff at the other, and second to include DH repulsions. Simulations of these librating chains are then compared with an analytical limit, called the planar-quadratic or pq approximation and based on the assumption that dihedral angle fluctuations are small. Simulations are also compared with the OSF result based on smooth curvature of elongated chains. The two analytic results, pq and OSF, both intended to describe stiff chains, are always different from each other, but either may agree with the simulations for appropriate limiting values of the parameters. The OSF result in particular is not quantitatively valid for all stiff chains, but only for the continuous, wormlike chain for which it was derived.1,2 The two analytic results for stiff chains and the simulation results for all degrees of flexibility concur in the asymptotic linearity of P in D2 but may differ greatly in the behavior of bond vector correlations and in the behavior of P versus D at small D. Flexible chains with different parameter values may show dramatic differences in such properties for small D, but all agree with the OSF power law. The differences can all be rationalized qualitatively on the basis of the metamodel of Khokhlov and Khachaturian21 (KK), but their quantitative results rely heavily on Gaussian models and the notion of stable blobs for detailed predictions, and these have little relevance to the present work.
10.1021/jp1005742 2010 American Chemical Society Published on Web 02/11/2010
3186
J. Phys. Chem. B, Vol. 114, No. 9, 2010
Fixman
In the remainder of the Introduction the models that have been studied by theory or simulation, rarely both, will be briefly surveyed and the rest of the paper will be outlined. Lengthier reviews of related material may be found elsewhere.4-6,17,13 Our main emphasis in the Introduction is on the models themselves and their suitability for simulations. Theoretical methods or approximations that have been used to study the models are mentioned incidentally. All of the models are regarded as physically sensible; none pretend to realistic depictions of actual chains. 1. OSF Model and Controversy. The debate has concerned one or the other of the following equations:1,2
P ) P0 + aD2
(1.1)
P ∼ P0 + aD2
(1.2)
The equality was derived for a thin, stiff wormlike chain,3 with P0 being the persistence length at vanishing Debye length or vanishing charge, and with a given by
aOSF ) q2B/4
(1.3)
where q is the number of elementary charges per unit length along the backbone contour of the chain and B is the Bjerrum length (about 7 Å for aqueous solutions at room temperature). The OSF approach was basically very simple. The electrostatic bending energy was calculated for a bent rod with uniform curvature and added to the intrinsic bending energy. The usual connection between bending energy and P in the theory of wormlike chains was then employed to derive eq 1.1. Any complication of the OSF model, such as finite thickness or consideration of discrete rather than continuous chains, will introduce a small length scale which D must greatly exceed before a simple relation could be expected between P and D, and then only if P is much less than the contour length. Thus D . 1, expressed in units of the bond length, is a minimum requirement for eq 1.1 to be applicable to a discrete model and the debate then concerns the validity of eq 1.2. Neither of the parameters a or P0 need retain its original interpretation for more complicated models, but the important issue is usually considered to be the power of D that occurs in the asymptotic eq 1.2. Increasing errors in the OSF calculation must be anticipated as the flexibility increases, and Barrat and Joanny have proposed4
P0 >
1 q2B
(1.4)
as a criterion for the validity of eq 1.1. Their derivation is explicitly limited by the condition that P . P0 and therefore q cannot be too small, and the applicability of this criterion to the present simulation results is uncertain because of the difference in models. 2. Continuous Chains of Gaussian Type. Most of the theoretical work has been based on a Hamiltonian that describes a pseudowormlike chain of coupled Gaussian bonds. This corresponds to the continuous version of the Harris-Hearst model,22 and it is therefore a model that is not well adapted to a realistic treatment of effects important at small length scales, such as those that determine the persistence lengths of flexible chains, or to simulations of discrete chains. Manghi and Netz13 have compared the different theories in this category through 2004, all of which are variational in nature, and concluded that
failure to recover the OSF power law was due to “the improper use of a cutoff at small length scales.” Manghi and Netz affirm the robustness of variational treatments, but this affirmation responded to a history of variational failures due, apparently, to delicate requirements that were not observed. The analytical results developed below may help to understand this delicacy. Comparisons between theoretical results for this model and simulations of discrete approximations to it have not been found. 3. Freely Jointed and Freely Rotating Chains. Freely jointed chains with fixed bond lengths5,6,14,23 and Gaussian bonds6,23 have been simulated, and the later work confirms the asymptotic OSF power law but there are no quantitative, theoretical predictions for the models. Polymer micelles modeled as freely rotating chains of charged hard spheres have been simulated24 and agreement with OSF theory was claimed, but in later work25 involving a direct study of angular correlations a somewhat lower power was reported. In neither paper was a value for the bond angle found, so the closeness of the model to a wormlike chain and the expectation of quantitative agreement with eq 1.1 cannot be determined. 4. Librating Chains. This is the model studied here. It has the same backbone structure as the freely rotating chain but adds terms Rφi2 and βφiφi(1 in the phase space of dihedral angles φi. The analytical result for P is subject to the conditions R > 0 and |β| < R/2, and agrees well with the simulation results only for large R, as would be expected from the derivation presented below. Simulations could in principle be executed for any values of R and β, and of course the freely rotating chain model is duplicated if R and β both vanish, but the only analytical result available for that limit is from a modified OSF calculation described below. 5. Bendable Rods. Simulations have been reported for freely jointed chains with fixed23 bond lengths or lengths constrained by a FENE potential,16 both models having an added potential controlling the bond angle and allowing adjustment of chain stiffness. The former work23 reported P linear in D for flexible chains and linear in D2 for stiff chains. The latter16 reported one P inferred from the decay of angular correlations on short length scales which is linear in D, and a second P describing the decay on long length scales linear in D2. This behavior will be compared below with the present findings. 6. More Complicated Models. Electrostatic models more sophisticated than those based on DH theory have already been cited,7-10 along with one example11 of the structural elaborations that might be needed for biopolymers such as DNA. These complications are doubtless relevant to the study of many real polyelectrolytes, but the object of this work is to clarify the theory of relatively simple models that are closely related to the OSF model. For this purpose the complicated models are merely a distraction and will not be considered further. The present work provides no reason to dispute the conclusion that the OSF calculation fails qualitatively for these more complicated models. 7. Metamodels. The blob picture of Khokhlov and Khachaturian21 (KK) is the closest thing we have to a metamodel that makes predictions for a class of detailed models similar to those considered here. KK hypothesize that the OSF result may be applied to a flexible polyelectrolyte chain that is treated as a wormlike chain of blobs, each a diffuse cloud of charge. The qualitative implications of this picture, that the minimum P is the value given by the OSF equality and that observations of larger P may be interpreted as due to blobs or curling in the space of angular correlations, are employed below. That POSF is a minimum P for a class of models is only a conjecture that might be provable; it is true for all results reported here.
Electrostatic Persistence Length
J. Phys. Chem. B, Vol. 114, No. 9, 2010 3187
The geometry of the chain model and the correspondingly revised OSF calculation are presented in section II. The potential energy calculation for arbitrary conformations is taken up in section III, the calculation of bond vector correlations and persistence length in section IV, the results and comparisons in section V, and the discussion in section VI, respectively. II. OSF for the Librating Chain The structural model consists of a discrete chain of bond vectors bi, each having unit magnitude and making angle θ with its neighbors; bi · bi+1 ) cos θ, so θ is the supplement to the usual bond angle. A dihedral angle φi measures the rotation of bi+1 and the following chain around the axis bi, and the sequence of all φi fully specifies the chain conformation. A positive quadratic form composed of terms Rφi2 and βφiφi(1 adjusts the stiffness and curling of the chain, and the minimum electrostatic energy of the charged chain is presumed to be obtained for the stretched out configuration in which all φi ) 0, that is, for the all trans conformation. In the all trans state the beads occupy the two edges of a ribbon, with the bond vectors crossing alternately from one edge to the other. For a consideration of deformations it is helpful to suppose that the stretched ribbon is oriented along the x-axis and lies flat in the xy plane. As θ increases from 0, the width of the ribbon increases from 0 and reaches a maximum width of unity as θ approaches π, since |bi| is chosen as the unit of length. The deformations of minimum electrostatic energy must be those for which all φi are near 0, and the most easily visualized deformations occur if all φi have the same small magnitude. If all φi are equal, e.g., φi ) φ, then, with the customary definitions of coordinate systems, the ribbon follows a straight line near the xz plane, rising slightly upward from the x-axis toward positive z, or falling downward toward negative z, depending on the sign of φ. If, instead, the φi have alternating signs with φi+1 ) -φi, then the ribbon curves in the xz plane and remains smooth except for a uniform wrinkling of small amplitude. This smooth, circular deformation of the ribbon within a plane, or actually within a thin disk, seems to be as close to the locally circular arcs of a rodlike chain as can be obtained for the different backbone geometry and gives an obvious modification to the coefficient a in eqs 1.1 and 1.3; that is, q must be reinterpreted as the linear charge density measured along the x-axis of the ribbon rather than along the zig-zagging backbone. Distances between charges measured along the ribbon axis x are decreased by a factor cos(θ/2). Assuming that each bead has λc1/2 elementary charges, then q ) λc1/2/cos(θ/2); a further division by cos(θ/2) is necessary to reflect our measurement of P in units of beads or bonds rather than in length along the ribbon, giving the modified eq 2.1 for the librating or freely rotating chain as
aOSF ) λc
B 4 cos3
θ 2
(2.1)
III. Energetics and Normal Modes A. Backbone Energy and Modes. Real chain polymers typically have several potential wells for each φi, among them the deepest well centered on the trans position. To have a simple potential that permits some analytical progress, a single quadratic well centered on the trans position has been used here, along with a quadratic, optional coupling between adjacent dihedral angles, that is,
N
UB )
N
∑ ∑ ub(|i - j|)φiφj
(3.1)
i)1 j)1
with the coefficients ub(|i - j|) taken to vanish if |i - j| > 1. The more general form of eq 3.1 will be used below for a quadratic approximation to the complete potential including electrostatic repulsion. It is of course a defect that eq 3.1 does not conform to the physical requirement for periodicity in each φi, but quadratic forms do permit analytical progress if the interval (-∞, +∞) is allowed for each dihedral angle. If the quadratic potential barrier is sufficiently high that values near |φi| ) π rarely occur, as they are for most of the parameter sets that have been studied, then little harm should be anticipated. The indexing of the chain is defined as follows. The beads are numbered from 0 to Nb, the bond vectors from 1 to Nb, and the dihedral angles from 1 to N, where N ≡ Nb - 1. There is no bond following bNb, so φNb is not needed. And φ1 serves only to rotate b2 around b1 relative to a plane defined by b1 and an axis vector in the external frame, so the magnitude of φ1 has no effect on the magnitudes of any relative bead distances rij; it is retained to simplify the indexing. The φi are measured in radians, the unit of energy is taken to be kBT, and all quantities in eq 3.1 are dimensionless. The N × N matrix of coefficients Uij ≡ ub(|i - j|) in eq 3.1 has a simple, tridiagonal form, the diagonal elements all being R ≡ ub(0), and the off-diagonal elements all being β ≡ ub(1). Therefore the matrix U with elements Uij is trivially related to the unit shift operator U0 of Morse and Feshbach26 by
U ) R1 + βU0
(3.2)
where 1 is the N × N unit matrix and U0 is the tridiagonal matrix with vanishing diagonal and unit off-diagonals. The orthonormal eigenvectors26 of both U and U0 are the columns (or rows) of the N × N orthonormal matrix E,
Eij ≡
( N +2 1 )
1/2
( Nijπ+ 1 )
sin
(3.3)
and the eigenvalues λbk of U for k ) 1, 2, ..., N are
( N kπ+ 1 )
λbk ) R + 2β cos
(3.4)
If β lies outside the interval (-R/2, +R/2), some of the λbk will go negative, at small k if β < -1/2R or large k if β > 1/2R; the corresponding librations will then be unstable and the quadratic form will not serve as an adequate physical model of the uncharged chain. Positive values of β favor negative correlations between adjacent dihedral angles and increase the flexibility of the chain for a given value of R; negative values of course do just the opposite. In alternative language, positive β helps the formation of blobs or crumpling and negative β helps their dissolution. Construction of normal modes qi with the definition
qi ≡
∑ Eijφj
(3.5)
VB ≡
∑ λbnqn2
(3.6)
gives
3188
J. Phys. Chem. B, Vol. 114, No. 9, 2010
Fixman
and the Jacobian of the transformation from angles to modes is unity, so the normal modes or dihedral angles serve equally well as the phase space of the chain, as long as the phase space of each is infinite or is so approximated, as is done here. The inverse of eq 3.5 is φi ) ΣEijqj since E is orthogonal. B. Electrostatics and the Quadratic Approximation. The potential energy VE due to electrostatic interactions is given in the Debye-Hu¨ckel approximation by
VE )
∑ uDH(rij)
(3.7)
energies associated with the modes qk, and second, they are not really needed for that purpose either, as will now be discussed. The preceding direct calculation of ue(|m - n|) is straightforward but was used only to confirm an alternative method that follows from eq 3.8 and preceding assumptions, and which was much faster in practice for long chains. The coefficients ue(|x|) are short-ranged functions of the index x, and are in fact taken to vanish for |x| > Nx, so the quadratic form VQ will be diagonalized by the normal modes defined in eq 3.5, given that Nx , N and that end effects are neglected. Specifically,
i m rotated as a rigid body by φm radians around the axis through beads m - 1 and m, that is, around the axis bm. VQm is the change in electrostatic energy due to the increase in φm and is easily calculated numerically, so ue(0) can be determined. If two dihedral angles, e.g., φm and φn with m * n, are both nonvanishing, then Q Vmn ) 2ue(|m - n|)φmφn + ue(0)(φm2 + φn2)
(3.10) and now the chain consists of three flat ribbons adjoined and Q is also easy to calculate twisted relative to each other. Vmn numerically, and along with it the coefficient ue(|m - n|) for all required values of |m - n|. The exponential decay of the coefficients with a decay constant proportional to κ|m - n| justifies a truncation of ue(|m - n|) at a maximum |m - n| ) Nx with an Nx value proportional to D, as is further discussed in the Simulations section. Equations 3.9 and 3.10 will be described as an analytical calculation of the quadratic expansion, but the formulas have never been written out for two reasons. First, the formulas would be very bulky and are needed only for numerical computation of the
∑ λekqk2
(3.11)
k)1
where Nx
λek )
∑
( Nkmπ + 1)
ue(|m|) cos
m)-Nx
(3.12)
Equation 3.4 is obviously a special case of eq 3.12. Since the Monte Carlo method requires routines that allow a calculation of VE for an arbitrary sequence of φ’s, one need only specify that a single qk of small amplitude is excited, calculate the φ’s for the prescribed qk, and then use those same numerical routines to make VE and λek. The resulting λek are smooth functions of k, as will be illustrated below, and it proved quite feasible to interpolate a relatively sparse subset with sufficient accuracy. In the quadratic approximation, although not in the more general Monte Carlo simulations, the qk are independently distributed Gaussian variables. With
λk ≡ λbk + λek
(3.13)
it can be readily confirmed that
λk )
〈φiφj〉Q )
1 2〈qk2〉Q 1 2
∑ k
EikEjk λk
(3.14)
(3.15)
where 〈 〉Q designates an average formed in the quadratic approximation in which VQ replaces VE. Monte Carlo simulations using VE give averages designated as 〈 〉M, which differ from the formally correct average 〈 〉 only through finite sampling errors. Comparisons between 〈φiφj〉Q and 〈φiφj〉M will be provided below. IV. p-Bonds and the Planar Approximation In this work the persistence length P is inferred from the correlations of pair bond vectors pi (p-bond vectors), defined by
pi ≡ bi + bi+1
(4.1)
Their normalized correlation is
Cn ≡ p-2〈pi · pi+n〉M
(4.2)
Electrostatic Persistence Length
J. Phys. Chem. B, Vol. 114, No. 9, 2010 3189
where p is the magnitude of any pi, namely p ) 2 cos(θ/2) in units of the bond length |bi|, and i and i + n are chosen far from the chain ends with n g 0. The pi are all colinear and parallel to the x axis when the chain is fully extended, and their correlations do not exhibit the oscillations in 〈bi · bi+n〉M that are due to fixed bond angles. An exponential decay of Cn,
Cn ∼ I exp[-n/P]
(4.3)
is expected for discrete chains only for large n, with a front factor I less than unity. Note that P here measures the decay of correlations with increasing n rather than increasing length along the straight ribbon edges, as discussed above for the modified OSF result. There is no need for further consideration of the p-bonds in the simulations. The construction of all b’s from the set of φ’s that define a chain configuration is well-known,27,28 and a previous summary29 of relevant formulas was used without change in the simulations. On the other hand, an analytic or numerical calculation of Cn based on the quadratic approximation is wanted for an accurate comparison with simulations, and for a better understanding of the approximation, and the methods developed for rotational-isomeric state models27-29 are not immediately useful for this purpose. An approximate expression for Cn is therefore motivated by the same considerations leading to the modified OSF result, eq 2.1, as follows. A more formal approach may be found in the Appendix. It gets no further but seems more likely to lead to quantitative estimates of errors and perhaps to a refined result. Equation 2.1 was derived on the assumption that the only relevant deformations of the discrete model are smooth bends into circular arcs, similar to those of a wormlike chain. The actual deformations are simpler in one respect: the wormlike chain in its elongated state can be bent in any direction but the discrete, librating chain cannot. In its elongated state the discrete chain is a flat ribbon in the xy plane, and no deformation lying entirely within that plane is possible for a chain with fixed bond lengths and angles. If all dihedral angles are small and if the chain is not too long, the rotated p-bonds will lie in the xz plane. A sufficiently long chain in this context is measured by the value of n sufficient for Cn to reach the exponential decay of eq 4.3. A plausible conjecture based on this picture is that the angle between the first and last p-bonds in the sequence is a linear combination Φn of all φj in the sequence, and an examination of the first few possibilities pi · pi+1and pi · pi+2 leads to the further conjecture n
Φn ) sin
∑
θ (-1)jφi+j 2 j)1
n>0
(4.4)
and therefore
Cpn ≡ 〈cos Φn〉M
Cpq n ≡ 〈cos Φn〉Q
is used, where 〈 〉Q indicates an analytic average taken with the quadratic approximation to the potential energy, in contrast to the exact or Monte Carlo average 〈 〉 or 〈 〉M implied by eq 4.5 for Cnp. Since Φn is a linear combination of Gaussian variables, eq 3.15 gives
1 2 2θ Cpq n ) exp - σn sin 4 2
(
)
(4.7)
where
Γnk2 λk k)1 N
σn2 ≡
∑
(4.8)
and i+n
Γnk ≡
∑
(-1)jEjk
(4.9)
j)i+1
The sampling of p-bonds was indicated in eq 4.2 and the following equations to start with some pi far from the chain ends. In practice the value of i was chosen to be (Nb - nx)/2, were nx is the largest value of n for which the correlations were computed. A few purely numerical properties of Cnpq inferred from eqs 4.7 and 4.8 are worth examining before simulation results are considered. First, the expectation that σn2 is proportional to n is fullfilled only if n is large, typically larger than 5-50, which may be regarded as the size of a “blob”; that size depends very much on the stiffness parameters (R, β). Second, the mode energy λk ≡ λbk + λek according to eq 3.13, is dominated by the backbone contribution λbk for much of the range of k, the electrostatic part λek increasing first gradually and then steeply only as k nears its maximum value of N, and the maximum λNe turns out to be proportional to D2. This behavior is illustrated in Figure 1 which, for reasons considered next, deals only with modes having wave numbers very close to the maximum value. The scaled values λ45(k) for D ) 45 decrease to about 10-3 for small k, with λ20(k) for D ) 20 being about 5 times larger; both curves are computed with λc ) 0.1 and θ ) 70°. Both curves use the parameters given for Run 8 of Table 1 but do not depend on simulation runs. It might be thought that the high energy modes at large k would make little contribution to σn2 or Cnpq, but that suspicion would be unfounded; it turns out at large k that Γnk2 is increasing much more rapidly with k than is λk. To illustrate this behavior, we consider the partial sums
(4.5)
Γnl2 λl l)1 k
fn(k) ≡ in the planar approximation. If n ) 0, the sequence is empty and Φ0 ) 0. The Appendix gives a formal derivation of eq 4.5 based on the assumption that all φ’s in the sequence are small. Comparison of the planar approximation with the original eq 4.2, both based on Monte Carlo averages, allows the planar approximation to be tested in isolation from the quadratic approximation to the electrostatic energy. If the quadratic approximation is now employed, then the notation
(4.6)
∑
(4.10)
with σn2 ) fn(N). For large n the ratio Fn(k) ≡ fn(k)/fn(N) stays near zero until k is very close to N and then rises very rapidly to unity according to Figure 1, where results are shown for D ) 45 and λc ) 0.1, again based on Run 8 of Table 1. It is is evident that the correlations between p-bond vectors separated by n bonds, and therefore also the persistence length inferred
3190
J. Phys. Chem. B, Vol. 114, No. 9, 2010
Fixman V. Simulations
Figure 1. In the two upper curves λD(k) ≡ λek/D2 vs k/N is shown for the values of the Debye length D given in the figure, the ordinate being scaled by a factor of 3. k is the mode index, and N is the total number of modes. In the two lower curves Fn(k) is the partial sum over modes that gives ln Cpq n , summed up to a maximum mode index value k and normalized to Fn(N) ) 1. λc ) 0.1 for all curves in this figure, and θ ) 70°. Parameters are those from Run 8 of Table 1.
from the correlations, are dominated by a progressively smaller fraction of the available modes the larger is n, and these are the modes of highest energy. It is therefore not surprising that a successful application of variational methods to the calculation of persistence lengths, such as that devised by Manghi and Netz,13 requires a careful exclusion of low energy modes from the variational treatment. Asymptotic properties are determined by relatively small parts of the available phase space.
A. Methods and Parameters. The simulations were based on the Metropolis30 Monte Carlo algorithm. Equilibration was initiated with the stretched out chain and attempted moves ran through all dihedral angles in sequence. The magnitude of moves was adjusted during equilibration to give an approximately 50% acceptance. Random numbers were obtained from the Numerical Recipes31 generator ran2. In addition to the physical parameters that define the system of interest, namely the number of bond vectors Nb, the bond angle supplement θ, the backbone potential parameters R and β, and the electrostatic parameters λc and D, a few additional parameters were introduced to speed up the calculation of the electrostatic energy given by eq 3.7. In principle, the double sum in that equation should extend over all distinct pairs of beads in the chain since all are assumed to be equally occupied by λc1/2 elementary charges, although bead pairs with |i - j| < 3 have relative distances fixed by the constrained bond lengths and angles and may be excluded. As was mentioned above, terms with |i - j| > Nx have been discarded, with Nx ) 10D for D > 32 and Nx ) 15D for smaller D, to speed up the calculations and incidentally to minimize excluded volume effects on the inference of persistence lengths. The total length of the chain was also truncated, typically at Nb ) 20Nx, although 16Nx was used for the larger D and Nb was never decreased below 2000, because a large interior region was wanted for the sampling of correlations. For |i - j| g 3 and j > i, the rij are functions of all dihedral angles from φi+2 through φj-1. Equations 3.1 and 3.7 then define the complete potential VP ) VB + VE as a function of all dihedral angles, and exp(-VP) is the weighting factor in the Monte Carlo simulations and gives averages designated 〈 〉M or just 〈 〉. Values of the algorithmic parameters seemed to be
TABLE 1: Table of Monte Carlo Simulations and Run Numbers Cited in the Text and Figuresa Run
θ (deg)
R
β
λc
00 PM
IM
aM
0 PM
P00 pq
apq
P0pq
aOSF
Dmx n
reps
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
70 70 70 70 70 70 70 70 70 70 70 70 70 15 70 70 70 70 70 70 30 15 15 5
128 64 16 2 0.25 128 64 16 2 64 16 2 0.25 2 8.2 8.2 8.2 8.2 8.2 2 2 2 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3.9 -3.9 -3.9 3.9 0 0 0 0 0
0.6 0.6 0.6 0.6 0.6 0.1 0.1 0.1 0.1 0.03 0.03 0.03 0.03 0.6 0.03 0.03 0.03 2.4 2.4 2.4 0.03 0.03 0.03 0.03
1617 765 193 24.7 3.09 1561 782 196 24.8 785 193 24.4 3.07 474 99.9 5.64 182 177 5.56 22.5 119 468 29 262
0.99 0.99 0.98 0.95 0.93 0.99 0.99 0.95 0.84 0.99 0.94 0.73 0.50 0.98 0.88 0.74 0.89 0.99 0.99 0.98 0.92 0.99 0.69 0.98
2.49 ( 0.03 2.28 ( 0.05 1.79 ( 0.02 1.56 ( 0.01 1.54 ( 0.02 0.44 ( 0.01 0.442 ( 0.006 0.384 ( 0.006 0.351 ( 0.005 0.150 ( 0.004 0.153 ( 0.002 0.176 ( 0.003 0.249 ( 0.005 0.87 ( 0.01 0.136 ( 0.002 0.158 ( 0.001 0.129 ( 0.002 6.00 ( 0.07 6.06 ( 0.13 5.8 ( 0.1 0.094 ( 0.001 0.057 ( 0.002 0.138 ( 0.002 0.060 ( 0.001
1650 ( 30 900 ( 60 400 ( 20 119 ( 7 120 ( 20 1600 ( 20 784 ( 10 230 ( 10 80 ( 10 778 ( 4 197 ( 2 53 ( 3 58 ( 5 507 ( 9 127 ( 3 100 ( 3 176 ( 3 400 ( 70 850 ( 160 300 ( 100 115 ( 1 469 ( 2 30 ( 2 259 ( 1
1556 778 195 24.3 3.04 1556 778 195 24.3 778 195 24.3 3.04 470 99.7 5.50 195 189 5.50 24.3 119 470
2.70 2.70 2.70 2.70 2.70 0.451 0.450 0.450 0.451 0.135 0.135 0.135 0.135 1.53 0.135 0.135 0.135 10.8 10.8 10.8 0.083 0.076
1550 771 187 18 -1.6 1555 777 193 23 778 194 24 2.8 466 99 4.7 194 167 -21 3 119 469
1.364 1.364 1.364 1.364 1.364 0.227 0.227 0.227 0.227 0.0682 0.0682 0.0682 0.0682 0.770 0.0682 0.0682 0.0682 5.46 5.46 5.46 0.0416 0.0385 0.0385 0.0376
456 456 456 456 456 455 528 455 525 459 459 394 455 457 606 605 606 457 455 456 396 398 398 397
1 1 2 1 1 4 1 3 2 2 3 4 2 2 2 4 1 1 1 1 2 2 1 4
a Except for the two columns on the right, the column headings are defined in the text. Dmx n is the largest value of D used to determine the slope aM from a plot of P versus D2, and the subscript n on Dmx n is the number of D values used in the plot. Those n values were chosen in order from the sequence 52, 45, 39, 32, 25, 20, 15, 10, 5, 2.5, 0, starting with the specified value of Dmx n . Thus each Run constitutes a set of runs with a sequence of D values, usually starting with D ) 0 and extending to Dmx n . The last column is the number of independent Monte Carlo repetitions, which differ only in the choice of starting seed for the random number generator. Probable errors are those reported by the least squares fits and are gross underestimates when the range of D values is small, as is further discussed in the text.
Electrostatic Persistence Length reasonable and other values were not extensively explored. Several values of θ were examined, namely 5°, 15°, 30°, and 70°. The largest angle is easiest for simulations, because it shows the greatest electrostatic effects for a given chain length, and is reasonably close to the proper value for carbon chains, so most of the results refer to 70° chains. The backbone potential parameter R was varied between 1/4 and 128. The smallest value was thought to be dangerously small, allowing easy access to |θ| > π, but no anomalies were actually apparent. The large upper value was found necessary to provide convincing evidence that the simulation results were indeed approaching the pq, analytical limit for chains with high charge densities λc. As previously remarked, values of β must be restricted to the interval (-R/2, +R/2) for librations to be stable in the absence of electrostatic interactions. Ideally, the expected exponential decay of p-bond correlations, as shown in eq 4.3, should have been confirmed for every set of parameters used in the simulations, but this was not considered necessary first because the dependence of ln Cn on n was always found to become linear at modest values of n, and second because the study of Cn at large n was often impossible. The largest chain length Nb used in the simulations was 9600, and significantly larger values were not practical with the available computational facilities. On the other hand, the value of P exceeded 104 for very large values of D and λc and end effects must be anticipated for large n. Even for short chains and small persistence lengths the study of n larger than a few times P becomes impractical because of statistical fluctuations. B. p-Bond Correlations. The dependence of Cn ≡ p-2〈pi · pi+n〉M on n is now examined for i in the interior of the chain, in comparison to the planar approximation Cnp and the quadratic, analytical approximation Cnpq. Cn from its definition and Cpn from the planar approximation eq 4.5 are calculated from Monte Carlo averages 〈 〉M based on an electrostatic energy which is complete except for the truncations discussed above. A third approximation Cnq may be constructed in which only quadratic terms are used in the energy but correlations are computed exactly from the φ’s, rather than invoking the planar approximation; this third approximation requires a separate but relatively rapid simulation and behaves, or misbehaves, similarly to Cnp in the few cases studied. The analytic pq approximation is seen in Figure 2 to benefit from compensating errors, but all errors are moderately small for the moderately large value R ) 16 used there. Worse agreement is found for the smaller R ) 2 used for Figure 3, and the discrepancy between Cn and Cpn becomes quite large as D increases. Irrespective of discrepancies between the analytic approximation and the simulations, the anticipated linearity of ln Cn in n is always recovered at moderate values of n and the persistence length P can be inferred from the exponential decay given in eq 4.3. Individual P’s obtained in this way were judged to have uncertainties in the range 2-5% on the basis of irregularities in the plots of major interest, namely P versus D2 or, more reliably, from different Monte Carlo runs on the same parameters. The following discussion will therefore be concerned only with the dependence of P on the various parameters, with probable errors in the slopes aM recorded in Table 1 from straight line fits. It should not be forgotten, however, that the other parameter I in eq 4.3, also tabulated in 0 for every Run, may drop appreciably Table 1 along with PM below its nominal value of unity for flexible and weakly charged chains. Those systems are generally found also to have appreciable curvature in plots of P versus D2 at small to
J. Phys. Chem. B, Vol. 114, No. 9, 2010 3191
Figure 2. Simulations and theory for Debye lengths of 5, 32, and 45, in units of the bond length. Lines are the analytic approximation ln p Cpq n , open circles are the planar approximation ln Cn, and filled diamonds are Monte Carlo results ln Cn. Angle θ ) 70°, charge λc ) 0.1, stiffness parameters R ) 16, β ) 0. Uses Run 8 in Table 1.
Figure 3. Simulations and theory for Debye lengths of 0, 15, and 39, in units of the bond length. Lines are the analytic approximation ln p Cpq n , circles are the planar approximation ln Cn, and diamonds are Monte Carlo results ln Cn. Angle θ ) 70°, charge λc ) 0.1, stiffness parameters R ) 2, β ) 0. Uses Run 9 in Table 1.
moderate values of D, and inference of the slope aM from a limited number of runs with large D may have much greater uncertainty than shows up in the least-squares fit. For the systems discussed in this work, Table 1 shows the 0 obtained from the simulation results for values of aM and PM 0 from the pq Cn at large n, the corresponding values apq and Ppq 00 approximation, and aOSF from the OSF approximation. PM is the simulation P for vanishing charge or D, i.e., in the absence 00 is the corresponding value of electrostatic interactions, and Ppq in the pq approximation. References to the Run number are to the first column of the table. It is seen that the agreement between aM and apq is made only slightly worse by the decrease in R from 16 to 2. Conversely, however, very large values of R are required to demonstrate that apq is approaching aM at large R. C. Stiff Chain Limit. We now consider the limiting behavior of the P versus D2 curves for increasing R and fixed β ) 0. In
3192
J. Phys. Chem. B, Vol. 114, No. 9, 2010
Figure 4. Persistence length P versus square of Debye length D, for backbone charge parameter λe ) 0.1 and for values of stiffness R equal to 128, 64, 16, and 2, as shown in the figure. Solid lines give the analytic approximation and the points are Monte Carlo results. The OSF slope is shown by the dot-dash line. Uses Runs 6, 7, 8, and 9 in Table 1, all having θ ) 70°.
Figure 5. Persistence length P versus square of Debye length D, for backbone charge parameter λe ) 0.6 and for values of R ) 128, 64, 16, and 0.25, as shown on the figure, and for the OSF calculation. The two solid lines at the top are computed in the pq approximation; the other two pq lines at the bottom are omitted but are parallel to those shown.The points are Monte Carlo results. The next two lines down, dotted, are simply linear fits to the Monte Carlo results, omitting the data for small D. The dot-dash line at the bottom is the OSF result. Uses Runs 1, 2, 3, and 5 in Table 1, all with θ ) 70°.
Figure 4 for λc ) 0.1, the Monte Carlo results are generally rather close to the quadratic approximation although systematic deviations are evident for small R and large D. These slowly disappear with increasing R. The intercepts IM shown in Table 1 seem to approach their nominal value of unity as R increases. Equation 1.2 is followed but the slopes apq are about twice aOSF, regardless of the value of R. If the charge density is increased to λc ) 0.6, the results displayed in Figure 5 indicate much larger deviations from the quadratic approximation, but again the deviations decrease with increasing R. Results for λc ) 2.4, a value much larger than is likely to be experimentally accessible, are shown in Runs 18-20 of Table 1 and indicate
Fixman
Figure 6. Persistence length P versus square of Debye length D, for stiffness parameters (R, β) ) (2, 0). Solid line Q gives the analytic approximation. Dotted line M is a least-squares fit to Monte Carlo results (solid diamonds) for the larger D. Angle θ ) 15°, charge parameter λc ) 0.6. The OSF result is shown by the dot-dash line. Uses Run 14 in Table 1.
that increasing charge density for fixed R depresses the slope aM toward aOSF, in agreement with the Barrat-Joanny criterion eq 1.4. In Figure 5 this approach of aM toward aOSF is seen even for the smaller charge density λc ) 0.6, but aM was never seen to decrease below aOSF. These observations seem qualitatively consistent with the KK treatment, but there are certainly no blobs that are sharply defined by the weak entropic and energetic forces considered here. Rather the electrostatic tension increases with increasing D, and slowly changes the initial slope of P versus D2. Only in the two limits of large R for fixed λc and large λc for fixed R, are the curves of P versus D2 seemingly straight lines for the full range of D, agreeing with the quadratic approximation in the former case and with the OSF result in the latter. Apart from the limiting cases considered above, an inference of asymptotic slopes is subject to systematic errors due to too small D that may be much larger than the probable errors shown in Table 1, which reflect primarily sampling errors in the simulations. At still smaller charge densities, some of which will be discussed below, the Monte Carlo results show P greater than given by the analytic approximation for small R, but decreasing slowly toward it as R is increased; the values of aM and apq for λc ) 0.03 reported for Runs 10-13 in Table 1 show that the ratio aM/apq falls from about 1.8 to 1.1 in this sequence of increasing R, from 0.25 to 64. Pushing to extremely stiff chains for small charge densities is difficult because the electrostatic contribution to the persistence becomes relatively small unless very large Debye lengths are used, and that requires longer chains and lengthier simulations than resources allowed. D. Bond Angle Dependence. 1. High Charge. The only significant differences between the rather flexible and highly charged chains shown in Figures 6 for θ ) 15° and 7 for θ ) 70°, both having (R, β) ) (2, 0) and λc ) 0.6, are the vertical shift upward and the decreased slope of Figure 6 relative to Figure 7. In both figures the persistence lengths lie closer to the OSF result than to the analytic approximation, as generally occurs at this moderately high charge and low stiffness. Both effects due to the decreased θ make an exploration of the stiff
Electrostatic Persistence Length
Figure 7. Persistence length P versus square of Debye length D, for stiffness parameters (R, β) ) (2, 0). Solid line Q gives the pq analytic approximation. Dotted line M is a least-squares fit to Monte Carlo results (solid diamonds) for the larger D. Angle θ ) 70°, charge parameter λc ) 0.6. The OSF result is shown by the dot-dash line. Uses Run 4 in Table 1.
Figure 8. Comparison of P versus D2 results for θ equal to 15°, 30°, and 70°, Runs 22, 21, and 12, respectively, all with (R, β) ) (2, 0). The 30° results are all displaced upward by 200. Filled diamonds are simulation results for P inferred from the correlation function Cn at large n. Solid lines are from the analytic pq approximation and dotdash lines give the OSF result.
chain limit more difficult and do not seem to have any compensating interest. Similar but smaller effects were found for 30° chains. A more interesting dependence on the magnitude of θ was found for low charge density. 2. Low Charge. For chains with a low charge density parameter λc ) 0.03, the results found for the dependence of P on D2 with θ ) 70°, 30°, and 15°, are shown in Figure 8. As usual in this work, P has been evaluated from the exponential decay of Cn at large n, but an attempt was made to infer a P characteristic of the decay of Cn at small n in view of recent simulations16 of stiff chains that examined the length scale dependence of P. In that work a linear dependence of P on D was found at short length scales and a quadratic dependence at
J. Phys. Chem. B, Vol. 114, No. 9, 2010 3193
Figure 9. P versus D2. All points are simulation results for charge λc ) 0.03 and stiffness R ) 8.2. Run 15 with β ) 0 gives the middle squares, and Run 16 with β ) 3.9 gives with simulations the bottom open diamonds and with the pq approximation the line Q2. The dashed line is a nonlinear least-squares fit to the open diamond points, P ) 1.53D1.48. Run 17 with β ) -3.9 gives the top filled diamonds and the pq approximation line Q1; the dotted line is a linear least-squares fit to the filled diamond points.
long length scales. In the present work neither the simulations nor the pq approximation allow such a simple characterization. Cn cannot be well fit as the sum of two exponents; rather the ln Cn vs n plots curve quite gradually to a straight line at large n. If some arbitrary maximum n, e.g., nmx, is used to for a straight line fit, the resulting P values have a quadratic dependence on D, although the values of P are of course smaller than those given from the asymptotic slopes. If nmx ) 1, then the inferred P values are very close to the value for an uncharged chain, although unrealistic values of λc give rise to as much as 20% increases of P with increasing D even for nmx ) 1. Very large θ values could probably do the same. The differences in backbone models, involving somewhat flexible bond lengths and angles in the other work,16 are presumed to be responsible for the different behavior. The asymptotic slopes in Figure 8, and the corresponding entries in Table 1, are qualitatively consistent with the conclusion that aM approaches aOSF as θ vanishes, and Runs 23 and 24 in Table 1 support this conclusion for the freely rotating chain. E. Effects of the β Stiffness Parameter. The effects of a change in β are very much dependent on the charge density λc; therefore, small and large charge density will be examined separately. 1. Low Charge. Figure 9 illustrates the effect of changing β for a moderately stiff chain with R ) 8.2 and a small charge density λc ) 0.03. The two nonvanishing values of β, the open diamonds at the bottom with β ) 3.9 and the filled diamonds at the top with β ) -3.9, are close to the extreme values permitted by the stability requirement |β| < R/2. Positive β at the bottom depresses the mode energies λk most at large k and negative β depresses λk mostly at small k. In physical terms, positive β favors alternating signs of succeeding dihedral angles and negative β favors constancy of sign, as will be illustrated below. The alternating signs for positive β correspond to increased curling of the chain and are similar in effect to the lowered stiffness caused by smaller R. The strong curvature in the P versus D2 plot for positive β can be easily fit with an empirical power law P ) 1.53D1.48, which might seem
3194
J. Phys. Chem. B, Vol. 114, No. 9, 2010
Fixman
Figure 10. P versus D2. All points are results for charge λc ) 2.4 and stiffness R ) 8.2. Run 18 with β ) -3.9 gives the filled diamonds from simulations, the dotted line as a least-squares fit to those points, and the pq approximation is labeled Q. Run 19 with β ) 3.9 gives the open diamonds from the simulations and the open circles from the pq approximation.
Figure 11. 〈φiφi+n〉 versus n for interior dihedral angles i, for Runs 16 and 17, all with R ) 8.2. The points are results from the Monte Carlo simulations and lines from the pq approximation. The filled diamonds are from Run 17 with β ) -3.9, λc ) 0.03, and D ) 45. The open diamonds are from Run 16 with β ) +3.9, λc ) 0.03, and D ) 45. The open circles are from Run 16 and D ) 0, which turns off electrostatic interactions and gives the same results as λc ) 0.
significant to some. Very likely to have much greater significance is the observation that the three sets of data covering almost the complete range of allowed β have similar slopes at large D, the obvious explanation according to the KK picture being that the increased tension in the chain due to increased D unravels the curls or blobs caused by positive β.32 Related effects occur in the 〈φiφi+n〉M and are discussed below. The set of five to six simulation points with largest D values may be fit to a straight line in D2 to excellent accuracy, as recorded in Table 1 and partially in Figure 9; the figure is too crowded to show least-squares fits to the open diamonds and squares at large D, but a straight edge overlay will easily confirm the claim. The strong disagreement between the open diamonds showing PM and the Q2 line showing Ppq in Figure 9, can be blamed almost entirely on failure of the quadratic rather than the planar approximation, because for this simulation Pp agrees very closely with the Monte Carlo values PM. 2. High Charge. The extreme charge density λc ) 2.4 in Figure 10 pushes the results close to the OSF line, and the change in β between the extremes of -3.9 and +3.9 has almost no visible effect. Large charge density alone is capable of decurling or deblobbing the polyelectrolyte. 3. Dihedral Angle Correlations. The results in Figure 11 pertain again to the low charge λc ) 0.03, as in Figure 9. Similar effects are found at high charge density but are much diminished in magnitude because large λc suppresses the oscillations in dihedral angles. Probably correlations over extended sequences of dihedral angles are important in the OSF limit. The smoothing out of oscillations in Figure 11 by electrostatic repulsions increases with increasing λc, and large charge densities give aM values close to aOSF and smooth curves of 〈φiφi+n〉M versus n. The smoothing of oscillations shown in Figure 11 for D ) 45 is essentially complete at D ) 32. The dramatic effects of deblobbing due to the change in β from β ) +3.9 to β ) -3.9 are seen in Figure 11 as well as Figure 9.
limit of increasing charge density q or λc for fixed θ and (R, β), or in the limit θ f 0 for any charge density and fixed (R, β). The planar-quadratic approximation pq developed above gives an apq about twice as large as aOSF, but the larger value is actually observed in the simulations only in the limit of increasing stiffness, through increased R, for fixed λc. On the other hand, the asymptotic eq 1.2 is confirmed for every choice of parameters that has been studied, within the error limits of our computational facilities. The simulations do confirm the implication that increasing q, all other parameters being held constant, drives P toward the OSF result, as is expected from the Barrat-Joanny criterion,4 but an increase in P0 may drive P either to OSF or to the pq result depending, respectively, on whether θ is decreased or R is increased. The planar-quadratic approximation has been quantitatively confirmed within its expected and limited domain of validity, namely under conditions that restrain all dihedral angles to small fluctuations from the trans conformation. The major control on the magnitude of these fluctuations in this work is the backbone potential parameter R, which must usually be increased to very large values to make the pq approximation quantitatively accurate. Under some conditions, deviations from the simulation results can be reduced and nearly quantitative agreement obtained by decreasing β to negative values, which reduces the crumpling otherwise entropically favored; see the example in Figure 9. On the other hand, a large charge density for small and moderate values of R always forces the simulation results for P to approach the OSF values; see Figures 7 and 10. Errors due to the quadratic expansion of VE are never conspicuously large for the 〈φiφi+n〉Q correlations, but these do not depend on the planar approximation and are calculated from eq 3.15, which reflects a broad sampling of mode properties across the spectrum; in other words, 〈φiφi+n〉Q gives only a weak test of the quadratic approximation and gives no test of the planar approximation. The inference of Ppq from eq 4.7 involves both the planar and quadratic approximations and is sensitive mainly to the very high wavenumber spectrum where the electrostatic contribution becomes dominant, unless the charge density and/or D become so small that the chain is little affected
VI. Discussion The continuous worm model and discrete approximations to it are quite similar but are not found to be identical: the value a ) aOSF in eq 1.1 is observed in the simulations only in the
Electrostatic Persistence Length
J. Phys. Chem. B, Vol. 114, No. 9, 2010 3195
by electrostatic interactions. The sensitivity of Ppq to a small part of the available phase space, and that part having high energy, makes variational calculations based on the minimization of free energy very delicate. The KK21 picture is qualitatively consistent with the results in that regions of large crumpling are associated with rapid increases in P above the OSF values, but the slopes of these curves always level off with increasing D and the curves may even become nearly parallel to those from the pq approximation; that is, aM may be only slightly larger than apq, although the simulations usually show a significant vertical displacement. This phenomenon is interpreted as a dissolution of blobs; likewise, oscillations in the 〈φiφi+n〉M versus n plots are smoothed out by the increasing D, and these effects seems certain to be part of the explanation of the common experimental inference that P is linear in Dx with a power x that is much less than 2 and may even approach 1. The main conclusions are, therefore, that the quantitative result for the wormlike chain found by OSF is by no means guaranteed even for similar models of stiff chains, but that the asymptotic linearity of P in D2 is found for every set of parameters that have been studied. The manner of approach to that limit and the magnitude of the coefficient of D2 may depend on the detailed parameters and modeling. Finally it should be mentioned that the values of D used here probably reach the largest values used in experiments. For a tetrahedral carbon chain with a bond length of 1.54 Å, the largest value D ) 60 used in the simulations corresponds to a Debye length of about 92 Å, corresponding to a univalent salt solution somewhat under 10-3 M.
(
cos γ sin γ 0 Tz(γ) ≡ -sin γ cos γ 0 0 0 1
Cn ) 〈eTx M(φi+1) M(φ3) · · · M(φi+n)ex〉
(A1)
where any φi measures a rotation around bi, not around pi. eTx is the row vector [1 0 0] and its transpose is the column vector ex. The 3 × 3 orthogonal matrix M(φ) is rather messy, as is construction of the direction cosines from the bk formulas. More simply, visualization of the defined transformations furnishes the formula
θ θ M(φ) ) Tx(π) Tz T (φ) Tz 2 x 2
( )
()
(A2)
Tx(φ) ) eφB
(
)
(A5)
where
(
0 0 0 B ≡ 0 0 -1 0 1 0
)
(A6)
This follows from expansion of the exponential and the observations that B2 ≡ -1yz, where
( )
0 0 0 1yz ≡ 0 1 0 0 0 1
(A7)
B3 ) -B, and so on. Gathering even and odd powers of B leads quickly back from eq A5 to eq A3. The required unitary transformation of Tx(φ) by Tz(γ) where γ ≡ θ/2 now gives
Tz(-γ)eφBTz(γ) ) exp(φBz)
(A8)
Bz ≡ Tz(-γ) BTz(γ)
(A9)
Bz ) B cos γ + C sin γ
(A10)
where
The matrix C is
(
0 0 1 C≡ 0 0 0 -1 0 0
)
(A11)
The remaining matrix Tx(π) occurs at the front of every M(φ) in eq A2, and in eq A1 Tx(π) will be replaced alternately, starting with the first M(φi+1), by either Tx(π) ) SS or the complex conjugate S*S*, where
( )
1 0 0 S) 0 i 0 0 0 i
(A12)
Since Sex ) ex the net effect of this rearrangement is that every factor exp(φBz) receives another unitary transformation, starting with [S exp(φi+1Bz) S*], [S* exp(φi+2Bz) S], and so on. The result in the exponentials is
where
1 0 0 Tx(φ) ≡ 0 cos φ -sin φ 0 sin φ cos φ
(A4)
and the product gives the same expression for M(φ) as does the more tedious calculation of direction cosines. Equation A2 has the additional advantage of suggesting the line of attack. The central matrix Tx(φ) in eq A2 can be rewritten as
Appendix A: Details of the Planar Approximation For specific calculations a reconstruction of the reference frames used for the bi is useful. The coordinate system associated with pi is specified by three unit vectors exi , eyi , and ezi defined as follows: exi points along pi, eyi points from the midpoint of pi to the end point of bi, and ezi completes a right handed orthonormal set. Then, just as for the bi,29 a sequence of n tranformations with 3 × 3 orthogonal matrices M(φ) composed of direction cosines between one frame and the preceding allow pi+n to be expressed in the frame of pi,
)
(A3)
SBzS* ) S(B cos γ + C sin γ)S*
(A13)
3196
J. Phys. Chem. B, Vol. 114, No. 9, 2010
Fixman
Cpn ) 〈cos Φn〉
or
SBzS* ) B cos γ - iF sin γ
(A14)
where
( )
0 0 1 F) 0 0 0 1 0 0
(A15)
S*BzS ) B cos γ + iF sin γ
(A16)
Likewise
The exponentials displayed in eqs A14 and the conjugate (A16) occur alternately in the product of M’s. To the first order of all φ’s that occur in the product of transformed M(φj), the nonvanishing commutator of B and F can be ignored in a Zassenhaus33,34 rearrangement of the exponentials, because it adds only higher powers of the φ’s. Then the coefficients of B can be gathered, e.g., on the right side of the product, and have no effect because Bkex ) 0 · ex if k˙ > 0. What remains is
Cn ) 〈eTx M(φi+1) M(φ3) · · · M(φi+n)ex〉
(A17)
Cn = 〈eTx exp[i(-φi+1 + φi+2...) sin γF]ex〉
(A18)
Cn ) 〈eTx exp(iΦnF)ex〉
(A19)
or
where n
Φn ) sin
∑
θ (-1)jφi+j 2 j)1
n>0
(A20)
One now observes that Fkex ) 0 · ex for odd k and Fkex ) ex for even k, so only even powers from the power series expansion of exp(iΦnF) in eq A19 contribute to the following result for the small φ or planar approximation, which is now designated Cnp:
(A21)
References and Notes (1) Odijk, T. J. Polym. Sci., Part B: Polym. Phys. 1977, 15, 477–483. (2) Skolnick, J.; Fixman, M. Macromolecules 1977, 10, 944–948. (3) In addition to assuming that the interactions have Debye-Hu¨ckel form, OSF assumed that the charges were uniformly distributed along a backbone of infinitesimal thickness, and that P . D; the last restriction justifies the treatment of electrostatic repulsions as a contribution to the bending energy of a wormlike chain. (4) Barrat, J. L.; Joanny, J. F. AdV. Chem. Phys. 1996, 94, 1–66. (5) Everaers, R.; Milchev, A.; Yamakov, V. Eur. Phys. J. E 2002, 8, 3–14. (6) Ullner, M. J. Phys. Chem. B 2003, 107, 8097–8110. (7) LeBret, M. J. Chem. Phys. 1982, 76, 6243–6255. (8) Fixman, M. J. Chem. Phys. 1982, 76, 6346–6353. (9) Netz, R. R.; Orland, H. Eur. Phys. J. E 2003, 11, 301–311. (10) Ariel, G.; Andelman, D. Phys. ReV. E 2003, 67, 011805-1-01180511. (11) Manning, G. S. Biophys. J. 2006, 91, 3607–3616. (12) Li, H.; Witten, T. A. Macromolecules 1995, 28, 5921–5927. (13) Manghi, M.; Netz, R. R. Eur. Phys. J. E 2004, 14, 67–77. (14) Nguyen, T. T.; Shklovskii, B. I. Phys. ReV. E 2002, 66, 021801. (15) Dobrynin, A. V.; Carrillo, J.-M. Y. J. Phys.: Condens. Matter 2009, 21, 424112-1–424112-11. (16) Gubarev, A.; Carrillo, J.-M. Y.; Dobrynin, A. V. Macromolecules 2009, 42, 5851–5860. (17) Dobrynin, A. V.; Rubinstein, M. Prog. Polym. Sci. 2005, 30, 1049– 1118. (18) Dobrynin, A. V. Macromolecules 2005, 38, 9304–9314. (19) Dobrynin, A. V. Macromolecules 2006, 39, 9519–9527. (20) Dobrynin, A. V. Curr. Opin. Colloid Interface Sci. 2008, 13, 376– 388. (21) Khokhlov, A. R.; Khachaturian, K. A. Polymer 1982, 23, 1742– 1750. (22) Harris, R. A.; Hearst, J. E. J. Chem. Phys. 1966, 44, 2595–2602. (23) Ullner, M.; Jo¨nsson, B.; Peterson, C.; Sommelius, O.; So¨derberg, B. J. Chem. Phys. 1997, 107, 1279–1287. (24) Cannavacciuolo, L.; Sommer, C.; Pedersen, J. S.; Schurtenberger, P. Phys. ReV. E 2000, 62, 5409–5419. (25) Cannavacciuolo, L.; Pedersen, J. S. J. Chem. Phys. 2002, 117, 8973– 8982. (26) Morse, P. M.; Feshbach, H. Methods of Theoretical Physics; McGraw-Hill: New York, 1953. See the material on the unit shift operator in chapter 1. (27) Eyring, H. Phys. ReV. 1932, 39, 746–748. (28) Flory, P. J. Statistical Mechanics of Chain Molecules; Interscience Publishers: New York, 1969. (29) Fixman, M. J. Chem. Phys. 1973, 58, 1559–1563. (30) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087–1092. (31) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical Recipes in Fortran, 2nd ed.; Cambridge University Press: Cambridge, U.K., 1992. (32) Comparison of asymptotic slopes in Table 1 for the runs in Figure 8 shows that aM > apq for the largest β. In the KK picture this means that deblobbing is not quite complete. (33) Wilcox, R. M. J. Math. Phys. 1967, 8, 962–982. (34) Scholz,D.;Weyrauch,M.J.Math.Phys.2006,47,033505-1-033505-7.
JP1005742
