Correlation between flexibility of chain-like polyelectrolyte and thermodynamic properties of its solution

Structural and thermodynamic properties of the model solution containing charged oligomers and the equivalent number of counterions were studied by means of the canonical Monte Carlo simulation technique. The oligomers are represented as (flexible) freely jointed chains or as a linear (rigid) array of charged hard spheres. In accordance with the primitive model of electrolyte solutions, the counterions are modeled as charged hard spheres and the solvent as dielectric continuum. Significant differences in the pair distribution functions, obtained for the rigid (rod-like) and flexible model are found but the differences in thermodynamic properties, such as, enthalpy of dilution and excess chemical potential, are less significant. The results are discussed in light of the experimental data an aqueous polyelectrolyte solutions. The simulations suggest that deviations from the fully extended (rod-like) conformation yield slightly stronger binding of counterions. On the other hand, the flexibility of polyions, even when coupled with the ion-size effects, cannot be blamed for qualitative differences between the theoretical results and experimental data for enthalpy of dilution.


Introduction
In more than 50 years of intensive polyelectrolyte research a large number of experimental and theoretical papers have contributed toward better understanding of these solutions (for review, see [1][2][3][4][5][6][7][8][9][10][11]). Nature and synthetic chemistry have provided polyelectrolytes of different shapes: they can be rod-like as, for example, DNA, flexible (chain-like) as are many of the synthetic polyelectrolytes. Their conformation in the aqueous solution depends on the chemical structure of the polyion, presence of low-molecular electrolyte, and polyelectrolyte concentration. Polyelectrolytes containing carboxylic groups, for example, poly(acrylic acid) or poly(methacrylic acid) can dramatically change their conformation in water; the polyion chain opens from a coiled to an extended configuration as the chain is neutralized. In addition, these solutions are characterized by strong intermolecular association [12].
Fully sulfonated poly(styrenesulfonic) acid and its salts represent a different class of polyelectrolytes. Due to the nature of sulfonic group, these (chain-like) polyions assume locally extended conformation in aqueous solution. Their properties are most often analyzed within the framework of the rod-like cell model [1][2][3] or Manning's condensation approach [13,14]. Both models are relatively simple, treating polyions as infinitely long, uniformly charged cylinders, or infinite lines. In addition to these simplifications, the system is treated as continuous dielectric, while interactions other than Coulomb are ignored. Despite severe shortcomings, these purely electrostatic theories yield important insights into the counterion-polyion interaction.
The basic parameter of classical polyelectrolyte theories is the linear charge density parameter λ, defined as (1.1) In equation (1.1), e 0 denotes an elementary charge, b is the length of the monomer unit carrying e 0 , r is relative permittivity of the system and as usual, k B T is a product of Boltzmann's constant and absolute temperature. There is a number of examples showing [3,10,[15][16][17][18] that in order to bring the Poisson-Boltzmann cell model calculation into agreement with experimental data for a given polyelectrolyte, λ eff should be introduced (λ eff > λ). For polystyrenesulfonic acid and its alkaline salts good agreement between the rod-like cell model calculations and the experimental data for the osmotic coefficient was obtained for λ eff ≈ 3.4, instead of the structural value [3], which is most often taken to be 2.80. This holds true for other physico-chemical properties [16,[18][19][20][21]. Although this fact is known quite some time, the explanation of its origins is unclear. One of the reasons often blamed for the disagreement described above is the fact that the polyions are treated as rigid cylinders, while realistic polyions are, depending on their molecular structure, more or less flexible. It is quite clear that a more advanced modeling of polyelectrolyte solutions should treat polyions as relatively flexible macromolecules that can adopt various conformations in solution. Seminal studies of short flexible polyelectrolytes with explicit counterions are 33603-2 due to Stevens and Kremer [22,23], who published molecular dynamics simulations of polyelectrolyte chains in solution. They used a freely jointed bead-chain polymer model to describe flexible polyions embedded in dielectric continuum; the sizes of monomers and counterions were given by repulsive Lennard-Jones potential. More recently other studies based on the similar model have been published [24][25][26]. This polyelectrolyte model explicitly handles long-range Coulomb interactions between the charges in solution and those fixed on the polymer backbone. In contrast to the older studies, which treat finite polyions as fully extended rigid cylinders [27,28], new works explore different conformations and correlate the charge "binding" and polymer conformations. Another important study is the work by Antypov and Holm [26]. These authors performed simulations of short rod-like polyions for the purpose of examining the most critical approximations of the rod-like cell model and the Poisson-Boltzmann theory.
Here we also have to mention an interesting recent contribution by Carrillo and Dobrynin [29], concerned with another class of polyelectrolytes in solution. These authors performed a study of hydrophobic (partially sulfonated) polyelectrolytes, simulating the transition between the extended and collapsed state [30,31]. This is an extreme case, which does not apply to highly charged and locally extended polyelectrolytes [3,[15][16][17][18][19][20][21], and is therefore not considered here. The solution properties of such partially sulfonated polyelectrolytes substantially differ from those we wish to model here.
In the present contribution we wish to address the effect of flexibility (departures from fully extended conformation) of a polyion chain on thermodynamic properties of solution. In order to get some quantitative information about the effect, we utilized the Monte Carlo method to explore two polyelectrolyte models, schematically shown in figure 1. According to the first one, i) the polyions were represented as freely jointed chains of charged hard spheres. The counterions were modeled as charged hard spheres and the solvent as a continuum with a dielectric constant of pure water under conditions of observation. The second model ii) differed from i) only by the description of the oligoions; the beads of each oligoion were ordered in a straight line and remained totally rigid during the simulation. Computer simulations were performed for three different degrees of polymerization; N m (number of beads) was 16, 32, and 48. In order to test the effect of the size of the counterion (hydration-dehydration effect), we also varied its radius r c . The two experimentally most thoroughly studied thermodynamic properties, i.e., the enthalpy of dilution (see, for example, [32,33] and the activity coefficient [9], were examined as functions of polyelectrolyte concentration. In addition, the pair correlation functions between the charged sites and counterions in solution, as well as the sites and ions themselves, were monitored and used to explain the calculated properties. The simulation results are discussed in view of the experimental results for these systems.

Models and simulations
In the present work, two different models were used. In the first model, flexible charged oligomers are made of tangentially bonded charged hard spheres of radius 2 Å, while in the second model the monomer units (beads) of the same size are ordered in a straight line and remain totally rigid (see figure 1 (b)) during the simulation. Counterions are in both models represented as charged hard spheres. Each monomer unit carries one negative elementary charge, while monovalent counterions are positively charged. The total number of counterions is thus equal to the total number of monomer units due to electroneutrality. We systematically varied i) the oligomer size, which was 16, 32, and 48 monomer units for both models, ii) the counterion's size, which was 1, 2, and 3 Å, and iii) the monomer concentration from 5 · 10 −5 to 0.2 mol/L.
Interaction between the charged sites is governed by Coulomb and hard sphere potentials: Since the beads, as well as the counterions, carry one elementary charge, Z i and Z j can be either +1 (for the counterion) or −1 (monomer unit). Note that the model used here is essentially the same as studied by Chang and Yehtiraj [24] and later by Bizjak et al. [25]. The model examined by

33603-3
Stevens and Kremer [22,23] is slightly different but the two models provide essentially the same insights.
The computer simulations were performed using the same code as before [25,34]. Simulations were carried out in a cubic box, whose size was bigger than the size of oligomer. Long range Coulomb potential was treated using the Ewald summation technique. The number of particles in the simulation box was from 1280 to 5760, depending on the chain lengths and concentration of monomers. For flexible polyions we used three types of moves. Chains were "moved" by using a) translation (the whole chain was translated for particular random vector), b) rotation of the shorter end of the chain and c) bend stretching. For rigid chains, the step a) and rotation of the chain around its center were utilized. The counterions were moved through the box using simple translation procedure. The equilibration run required 1 · 10 5 , and to obtain good statistics for the production run, 1 · 10 6 moves per particle were needed.

Results and discussion
Numerical results in this section apply to the concentration of monomer units c m = 0.001 monomol/dm 3 . Further, T = 298.15 K, where L B is equal to = 7.14 Å, which yields the linear charge density λ = 1.785. In this way, λ > 1, this is considered to be a moderately strong polyelectrolyte for which the Manning condensation [13,14] should take place. To orient the reader, λ for poly(styrene)sulfonic acid is 2.8, and for the DNA it is 4.2, while for another class of frequently studied polyelectrolytes, i.e. x, y-ionenes [18][19][20][21], it is from 0.673 to 1.437. Only fully charged oligomers -where each bead carries one negative elementary charge -were studied here. The effect of the polyion length was explored with varying number of segments N m from 16 to 48.

Pair distribution functions
Here we are interested in how an increase in the flexibility of polyions and the size of counterions affect the spatial distributions of particles in such systems. In the binary systems like these ones we have three pair distribution functions (pdf): the monomer-monomer g mm (r), monomer-counterion g mc (r), and counterion-counterion g cc (r) function. These correlation functions have been calculated for both models and for two (or in some cases three) values of the counterions radius, r c . Only the most interesting pair distribution functions are presented here.

Flexible oligoions
In figure 2 (a)-(c) we present the results for the intermolecular part of the monomer-monomer distribution function for N m = 16, 32, and 48, all for two different values r c (1 and 3 Å), as indicated in the caption. This result shows that smaller counterions (1 Å), due to their stronger electrostatic screening of polyions, cause the monomer units belonging to different oligoions to distribute at distances (continuous curve) smaller than larger ions 3 Å (dashed curve). The full monomer-monomer distribution function, which also contains the intramolecular part, is shown in figure 3 The monomer-counterion pair distribution functions, shown in figure 4 (a)-(b), reflect strong counterion-oligoion attraction. These results apply to N m = 48. The counterions having smaller radius are capable of closer approaching the oligoion, and thus the peak of the pair distribution function reaches higher values. The minimum of g mc (r) at large distances, i.e. at r ≈ 200 Å is an interesting new feature still unknown within the Poisson-Boltzmann cell model calculations. The minimum was also found for N m = 16 (and 32) but it was shallower and located closer to the origin.
The third pair distribution function, shown in figure 5 (a)-(c), describes the counterion-counterion correlation, g cc (r). In several recent papers [35][36][37][38], this distribution function has been extensively studied for the cylindrical model with infinitely long polyion. A strong increase of concentration of equally charged counterions in the vicinity of the polyion forms a basis of the catalytic effect [39,40]. The shape reflects the fact that the highest probability of finding two counterions 33603-4  is just on the opposite sides of the charged bead. Such distributions of counterions are not limited to linear polyelectrolyte and may also appear in highly charged micellar solutions [41].

Rigid versus flexible oligoions
First studies of short cylindrical polyions, with the purpose of exploring the effect of an increased length of the oligoion on the osmotic coefficient, were published long ago [27,28]. The Monte Carlo method and the Poisson-Boltzmann theory were used for this purpose. These calculations were performed in the cell model approximation and the first simulation of an isotropic solution of 33603-5

33603-6
rigid-chain polyelectrolyte seems to be due to Antypov and Holm [26]. These authors, however, did not present any spatial distribution functions; the pair correlation functions for our system are shown next. In particular, we are interested in the comparison of the correlation functions obtained for flexible and rigid oligoions, which is presented next.
First in figure 6 (a)-(c) we display the result for the interparticle part of the monomer-monomer pair distribution function as obtained by the flexible and rigid oligoion models. The results show that monomers of flexible oligoions are generally distributed at distances larger than those of rigid oligoions. The other distribution functions are shown in figures 7 and 8. First, the monomer- counterion distributions are displayed in figure 7. As we see the distributions differ qualitatively; the descent of the function obtained for rigid polyions is steeper at the beginning and forms a much shallower minimum. At a certain distance, the two functions intersect each other. Up to this distance, the concentration of counterions near the rigid polyion is significantly higher. Different shapes of the pair distributions, of course, affect the fraction of counterions calculated within a certain distance from the polyion.
In figure 8 we show the counterion-counterion, g cc (r), distribution functions for both cases. Similarly, we observe significant differences between the functions obtained for rigid and flexible polyions. The correlation function for rigid polyions (dashed line) runs below the one for flexible polyions (continuous line), indicating weaker correlation between equally charged ions in the vicinity of the polyion. This observation is only true to a certain distance where, as shown in figure 7, the two distribution functions intersect. As found in previous studies of infinitely long and rigid polyions [36], smaller ionic size yields stronger correlation between the counterions in the electrical double-layer. This counterion-counterion correlation is neglected in the Poisson-Boltzmann theory but as shown before [35] (see also the references therein), this is not considered to be a serious source of error when monovalent counterions are present in solution. For divalent counterions this is not true any more. In such cases, the counterion-counterion correlation may be strong enough to cause 33603-7  clustering of the polyions; this effect cannot be accounted for by the mean-field Poisson-Boltzmann theory.

Conformational averages
The quantities, reflecting the average conformation of the oligoion in solution, are the end-toend distance and the average radius of gyration. Both properties have been for flexible oligoions studied before [25], but only for one size of the counterions. Here, the study is extended to the counterions with r c = 1 and 3 Å in order to quantify the effects of the counterion size on these averages.
The average values of the end-to-end distance, R ee (in Å) and the radius of gyration R g (in Å) are, as a function of the counterion size and for the degrees of polymerization N m equal to 16, 32, and 48, collected   An increase of the degree of polymerization increases R ee and R g . Also, larger counterions do not screen the repulsion between the monomer units (beads) as well as the smaller ones; the effect increases the end-to-end distance. It is of interest to compare the end-to-end distance for a certain oligoion, with its contour (fully extended) length. For N m = 48, the fully extended oligoion measures 192.0 Å, while the average end-to-end distance is only around 118.5 Å.

Enthalpy of dilution
For polyelectrolyte solutions the excess internal energy is not that important per se, however, the difference of U ex between two concentrations, can be approximated with the enthalpy of dilution, ∆ d H. This quantity can be accurately measured [17,32,33,[42][43][44][45], providing important insights into the nature of interaction in polyelectrolyte solutions. The measurements performed for various polyelectrolyte systems indicated important deviations from theoretical predictions (for review see, [44]). For example, the classical electrostatic theories predict exothermic effects upon dilution [3], while experimental data show that ∆ d H can be negative or positive depending on the nature of the counterion [33,45].
The results for ∆ d H ex are presented in figure 9 for different N m as a function of the − log 10 c m . The results for both flexible and totally rigid model of oligoions are presented in the same figure.
The most important observations are as follows: i) the enthalpies are, as found in our previous study [25], negative and the slope depends on the degree of polymerization; more heat is released upon dilution per mole of monomer units for shorter polyions, ii) the ∆ d H values are less exothermic for flexible polyions than for the rigid ones, iii) the effects of counterion radius are noticeable for solution of rigid polyions but much less pronounced in case of flexible polyions. In all cases the differences are small and, as our experience shows, they are within the uncertainties of experimental data. We can compare these results with those presented, for example, in figure 5 by Lipar et al. [17]. These authors have measured the heats of dilution of various alkaline salts of polyanetholesulphonic acid at 298 K. As we see both the trends and magnitude of the effect are correctly predicted by our models.

Activity coefficient of counterions
A comprehensive examination of the counterion activity in polyelectrolyte solutions was published by Wandrey et al. [9]. The study revealed the effect of the degree of polymerization on the counterion activity, the factor which was not evident from the theories treating polyions as infinitely long and rigid cylinders. Here we explore the effect of another factor of this kind, i.e., the polyion flexibility. Excess chemical potential (activity coefficient) of counterions µ ex has been studied in our previous work [25]. The novelty of the present work lies in the simulations that we performed for both i) flexible and ii) rigid models as well as in examining the effect of the size of the counterions.
As before, we use Widom's method to determine this quantity. In equation (3.1), ∆U represents the sum of pairwise potentials of the inserted fictitious counterion with all others counterions and monomer units (see, for example, [46]) Although the individual activity coefficient of counterion is, strictly speaking, not a well-defined thermodynamic property, still it is of theoretical interest. For the salt-free solutions in the cylindrical cell model approximation, this quantity is equal to the osmotic coefficient [47]. Earlier experimental determinations of counterion activity coefficient are due to Nagasawa and Kagawa [48], Oman and Dolar [3,49], and Joshi and Kwak [50]. More recent experimental results were, for solutions of aliphatic ionenes with variable degree of charging, obtained by Nagaya et al. [18] and by Minakata and coworkers [51]. From theoretical viewpoint, we have to mention the work of Manning [13]; his theory describes the individual ionic activities of electrolyte-polyelectrolyte mixtures in semi-quantitative agreement with experimental data. Recently Manning and Mohanty extended the theory to solutions of charged oligomers [52], emphasizing the the importance of the ratio of the Debye screening length L D to the contour length of the polyion.
The results of our Monte Carlo simulation are shown in figure 10. The reduced excess chemical potential of counterions (logarithm of the individual activity coefficient) exhibits quite a significant r c dependence. The counterions of smaller size yield more negative values of the excess chemical 33603-10

Conclusions
New simulations and theoretical results were presented for two models of polyelectrolyte solution. In one case, the polyions were pictured as a rigid one-dimensional array of charged hard spheres being totally rigid, while in the other example, they were represented as a freely jointed (flexible) chain. The results shown here strengthen and complement the conclusions arrived at in some previous papers, most notably in [25,26].
Monte Carlo simulation was utilized to calculate the pair distribution functions between various sites. The results for the monomer-monomer, monomer -counterions, and counterion-counterions distribution functions were obtained for the both oligoion models, varying the radius of the counterions from 1.0 Å to 3.0 Å. The site-site pair distribution functions indicate significant differences between the rigid and flexible model of the oligoion. In particular, the counterion-monomer pair distribution and monomer-monomer distribution functions can differ quite substantially if the counterions present in solution are small.
Among thermodynamic properties we choose to examine the enthalpy of dilution and the activity coefficient of counterions. These two properties, taking into account the correlation between the activity and osmotic coefficient, are widely used to characterize polyelectrolytes in solutions. Experimental results for the enthalpy of dilution very clearly reveal the ion-specific effects (see [33] and the references therein). Aqueous solutions of 3, 3-ionenes (cationic polyelectrolytes having quaternary ammonium groups on the backbone) with fluoride (strong kosmotrope) counterions yield an exothermic effect upon dilution in good agreement with theory and the simulation results presented here. On the other hand, the solutions with chloride, bromide, and iodide counterions consume the heat upon dilution (endothermic effect) in strong disagreement with the theory and simulations based on the continuum solvent models. The 3, 3-ionenes, studied in several papers [16,[19][20][21]33], have the charge density parameter λ = 1.437, which is not far from the value of the present model, with λ equal to 1.785. It is shown that the parameters used to describe the polyelectrolyte models studied here are not unrealistic. It is quite clear from the results presented in this paper (cf. figure 9) that neither the effect of flexibility of polyions, nor the size of the counterions (both in reasonable limits) can explain qualitative differences between the continuum-solvent theories and experimental results. As discussed elsewhere [33], the ion-specific effects discussed above, are closely  connected with the hydration properties of the counterions and charged group on the polyion.
The counterion activity coefficients in polyelectrolytes were studied in many papers and the results were thoroughly analyzed by Wandrey and coworkers [9]. The experimental data confronted with approximate theories. The analysis revealed important effects of the polyelectrolyte concentration and the chain length of the polyion. Measurements indicate that the activity coefficient of counterions decreases with the increasing degree of polymerization. This finding is consistent with our calculations. Our simulations, in agreement with measurements, predict a substantial ionspecific effect; the smaller is the ion the stronger is the ion binding [16,18]. In view of theoretical nature of our study, where we can control the model parameters, we have been able to examine the effect of the flexibility of oligoion on the counterion activity coefficient. Such effects cannot be easily studied experimentally. Our simulations show that flexible oligoions (or polyions) attract the counterions more strongly than the rigid ones. This finding supports the idea that λ eff > λ, should be used to successfully fit the experimental data.
In conclusion, the corrections due to the deviations from fully extended polyion conformations are generally not large. This is certainly true for highly charged polyelectrolytes examined here. The situation, however, may be different in case of weakly charged polyelectrolytes, or in the presence of large hydrophobic domains on the polyion. In such situations, transitions from the extended shape to the globular one are possible.