Langmuir 1999, 15, 4123-4128
4123
Monte Carlo Simulations of Ion Activities in Rodlike Polyelectrolyte Solutions† Takuhiro Nishio* and Akira Minakata Department of Physics, Hamamatsu University School of Medicine, 3600 Handa-cho, Hamamatsu 431-3192, Japan Received September 8, 1998. In Final Form: November 30, 1998 Canonical Monte Carlo simulations are applied to the investigation of the small ion distributions around a rodlike polyelectrolyte in aqueous solutions with or without added 1:1 salt. The mean and single ion activity coefficients are estimated and are compared with the experiments and the theories. The Monte Carlo cell systems are examined to measure the ion distribution adopting simple molecular models of the polyions, and the activity coefficients are obtained from the local ion concentrations at the cell boundary. The results are compared with the solutions of the Poisson-Boltzmann equation. On comparison with the experimental data in the literature, considerable agreement is observed in several cases. These comparisons provide the fundamental knowledge of our simple systems as well as of the essential properties of the polyion-small ion interactions. The applicability and the limitation of the method are discussed.
Introduction
Method
The colligative properties in the polyelectrolyte solutions at finite concentration are extensively treated with the cell model of the cylindrical coordinate.1,2 These properties are successfully interpreted by the limiting laws based on the counterion condensation hypothesis.3 However, several essential problems remain unsolved.4 Furthermore, the quantitative examination may not be carried out sufficiently even in the rodlike polyelectrolyte at finite concentration. The solution of the Poisson-Boltzmann equation is an effective way to deal with the cell system in an ideal case.5 The canonical Monte Carlo simulation method can be also used to investigate the mobile ion distribution and the thermodynamic properties. The advantage of the Monte Carlo method is the inclusion of nonideality in the real solution. The simulations of the potentiometric titration were successfully achieved for the (carboxymethyl)cellulose solutions using an isolated cylindrical cell system with a simplified molecular model.6 It seems that examination of the ion activity is valuable in applying a similar method. In the present study, ion distribution around a rodlike polyion molecule is evaluated with or without added monomonovalent salt by the canonical Monte Carlo method. The mean and single ion activity coefficients are estimated from the local ion concentrations at the cell boundary. The results of the simulation are compared with the experimental data in the literature and the solutions of the classical Poisson-Boltzmann equation. Two distinctive cell systems are examined for the Monte Carlo simulations. In the Poisson-Boltzmann calculations, the polyion charge is assumed to be smeared uniformly on the surface of the polyion rod. The applicability and the limitation of the method are discussed.
The mean ion activity coefficient in the PoissonBoltzmann cell model can be formulated in the case of mono-monovalent salt7
† Presented at Polyelectrolytes ‘98, Inuyama, Japan, May 31June 3, 1998.
(1) Katchalsky, A.; Alexandrowicz, Z.; Kedem, O. In Chemical Physics of Ionic Solutions; Conway, B. E., Barradas, R. G., Eds.; John Wiley & Sons: New York, 1966; p 295. (2) Katchalsky, A. Pure Appl. Chem. 1971, 26, 327. (3) Manning, G. S. J. Chem. Phys. 1969, 51, 924. (4) Ise, N. Adv. Polym. Sci. 1971, 7, 536. (5) Anderson, C. F.; Record, M. T., Jr. Annu. Rev. Biophys. Biophys. Chem. 1990, 19, 423. (6) Nishio, T. Biophys. Chem. 1996, 57, 261.
γ( ) [C+(R)C-(R)/CS(CP + CS)]1/2
(1)
where C+(R) and C-(R) are the local concentrations of counterions and co-ions, respectively, at the outer boundary of the cylindrical cell whose radius is R. CP and CS are the mean concentrations of the polyion and salt, respectively, in the whole cell volume, given by
CP )
NP 1 NAV πR2H
(2)
CS )
NS 1 NAV πR2H
(3)
where NP is the number of charged groups on the polyion, NS is the number of ions from excess salt, H is the height of the cell, and NAV is Avogadro’s constant. As a consequence of natural assumption, the single ion activity coefficients of counterion γ+ and co-ion γ- are expressed, as follows:8
γ+ ) C+(R)/(CP + CS)
(4)
γ- ) C-(R)/CS
(5)
The above definitions of the ion activity coefficients do not include the contribution of the interactions between small ions. Thus, the calculated value must be compared with the following “corrected activity coefficient” γc from the experiment
γc ) γex/γ°
(6)
where γex is the experimental activity coefficient and γ° is the activity coefficient of added salt.9,10 These relations (7) Marcus, R. A. J. Chem. Phys. 1955, 23, 1057. (8) Gue´ron, M.; Weisbuch, G. J. Phys. Chem. 1979, 83, 1991. (9) Wells, J. D. Biopolymers 1973, 12, 223. (10) Wells, J. D. Proc. R. Soc. London, Ser. B 1973, 183, 399.
10.1021/la981188r CCC: $18.00 © 1999 American Chemical Society Published on Web 03/13/1999
4124 Langmuir, Vol. 15, No. 12, 1999
Nishio and Minakata
Figure 1. Schematic views of two cell systems for the present Monte Carlo simulations: a, isolated cylindrical cell system (MC-c); b, hexagonal cell system (MC-h).
are assumed to be applicable to the canonical Monte Carlo simulation of the ion distribution. The comparisons are expected to provide the knowledge of the ion activity in the polyelectrolyte solutions. Rigorously speaking, ion activity coefficient should be determined as follows
ln γi )
(
)
∂(Fex/VkBT) ∂ni
(7)
T,V,n
j*i
where Fex is the excess free energy in the system, V is the volume of the system, kB is the Boltzmann constant, T is the absolute temperature, and ni is the number concentration of ith ion. Under the electroneutrality condition, mean ion activity coefficient is derived from this equation. The grand canonical Monte Carlo method is also useful to evaluate the mean ion activity directly.11,12 However, more complicated potential functions must be necessary for the detailed comparison of the activity coefficient. Furthermore, single ion activities cannot be estimated by these methods. The above-mentioned canonical Monte Carlo method seems suitable enough to extract the characteristics of the polyelectrolyte solutions for the present purposes with some reservations. A hexagonal cell system as well as a cylindrical cell system is attempted as a Monte Carlo cell system. The ordinary isolated cylindrical cell system (“MC-c” refers to the system) was previously applied to the simulations of potentiometric titration.13 This is made up of an objective cell with two external cells on both sides, as shown in Figure 1a.14 In the isolated cell system, however, inhomogeneity in the ion distribution appears in the region near the cell surface.15 The present hexagonal cell system (“MC-h”) is constituted by the packed 21 hexagonal cells (one central cell and the extended image over 20 surrounding cells), as illustrated in Figure 1b. The cell wall is impenetrable to the center of mobile ions. The interaction energy in the central cell is calculated taking account of (11) Mills, P.; Anderson, C. F.; Record, M. T., Jr. J. Phys. Chem. 1986, 90, 6541. (12) Jayaram, B.; Beveridge, D. L. J. Phys. Chem. 1991, 95, 2506. (13) Nishio, T. Biophys. Chem. 1994, 49, 201. (14) Mills, P.; Anderson, C. F.; Record, M. T., Jr. J. Phys. Chem. 1985, 89, 3984. (15) Wennerstro¨m, H.; Jo¨nsson, B.; Linse, P. J. Chem. Phys. 1982, 76, 4665.
Figure 2. Molecular models of (carboxymethyl)cellulose (CMC) of DS ) 1.0 (a) and pectate (DS ) 1.0) (b). Top view and the side view of one repeating unit are presented for each molecule with the volume from which the centers of the mobile ions are excluded.
all ions in the image cells. Although the long-range interactions cannot be included sufficiently, the surface effect can be almost eliminated in this multiple cell system. Similar systems are well studied for the osmotic pressure of DNA solutions by other researchers.16,17 Sodium (carboxymethyl)cellulose (Na-CMC) and sodium pectate are chosen as rodlike polyelectrolytes for the present study. The molecular model of CMC is constructed by an elliptic cylindrical backbone rod and spherical charged groups, as shown in Figure 2a.6 The number of the groups attached to a glucose unit of CMC is determined to be the same degree of substitution (DS) to the experiment. For the simulation of pectate, a molecular model of 31 symmetry is built up with a rod and spheres according to the X-ray diffraction analysis (Figure 2b).18 The degree of substitution of the sample is given to be 1.00 in the calculations instead of the experimental value 0.99, for convenience. The lengths of the monomer unit are defined to be 5.18 and 4.45 Å for CMC and pectate, respectively. The radius of carboxylic acid groups and small ions is fixed to be 2.0 Å. The Monte Carlo trial of the ion movement is carried out following the standard Monte Carlo sampling algorithm.19 The energy between two ionic charges is defined as the Coulombic interaction with a hard core repulsion. The small ion distribution is evaluated by summing up the ions in the coaxial cylindrical shell counters during a sufficiently long trial loop. The ion activity coefficient is determined from the resultant concentration in an outermost shell. The iterations are carried out until the mean error is below 0.3% in most cases. The solution of the Poisson-Boltzmann equation in the cylindrical cell system was obtained analytically in the case of no added salt.20,21 In the presence of excess salt, accurate numerical solution is achieved by the iterative (16) Nilsson, L. G.; Guldbrand, L.; Nordenskio¨ld, L. Mol. Phys. 1991, 72, 177. (17) Lyubartsev, A. P.; Nordenskio¨ld, L. J. Phys. Chem. 1995, 99, 10373. (18) Walkinshaw, M. D.; Arnott, S. J. Mol. Biol. 1981, 153, 1055. (19) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (20) Alfrey, T., Jr.; Berg, P. W.; Morawetz, H. J. Polym. Sci. 1951, 7, 543. (21) Fuoss, R. M.; Katchalsky, A.; Lifson, S. Proc. Natl. Acad. Sci. U.S.A. 1951, 37, 579.
Ion Activities in Rodlike Polyelectrolyte Solutions
Langmuir, Vol. 15, No. 12, 1999 4125
Figure 3. Counterion activity coefficients in salt-free NaCMC solutions. Dependence of Na+ activity coefficient upon polyion concentration without added salt. The DS value of CMC sample was 0.736 (ξ ) 1.014). Circle, experimental data; square, MC-c; +, MC-h results. Curves denote the Poisson-Boltzmann solutions. Polyion radii (a) are 4.5, 3.5, 2.5, 1.5, and 0.5 Å, from top to bottom.
Figure 5. Comparison with the data of Na-CMC-NaCl by Joshi and Kwak. Corrected activity coefficients in Na-CMCNaCl solution: γ-c, γ(c (closed symbols), and γ+c from top to bottom. Small symbols, experiments; large symbols, MC-h; lines, Poisson-Boltzmann (a ) 3.5 Å) results. Solid lines and circles, CP ) 0.01 N; broken lines and triangles, CP ) 0.025 N. The DS value of CMC sample was equal to 0.86 (ξ ) 1.185).
parameter, ξ, is given as
ξ ) lB/b
(8)
where lB is the Bjerrum length and b is the mean axial distance per group on the polyion. The Bjerrum length, which is taken as 7.135 Å, is defined as follows
lB ) e2/(4π0DkBT)
Figure 4. Comparison with the data of Na-CMC-NaCl by Rinaudo and Milas. Corrected activity coefficients in Na-CMCNaCl solution: γ-c, γ(c (closed symbols), and γ+c from top to bottom, depending on X1/2. Small symbols, experiments; large symbols, MC-h; lines, Poisson-Boltzmann (a ) 3.5 Å) results. Solid lines and circles for DS ) 1 (ξ ) 1.377); broken lines and squares for DS ) 2.5 (ξ ) 3.444). CP ) 0.002 N.
solution of the initial value problem for the onedimensional differential equation using Gear’s method.22 The activity coefficients can be calculated from the concentration at the cell surface. In the system, the polyion charge is assumed to be distributed uniformly on the surface of the cylindrical polyion rod. The solution of the Poisson-Boltzmann equation depends on the polyion radius, a, even in the same charge density. The temperature applied in the calculations is 25 °C, the same as in the experiments. The linear charge density (22) Byrne, G. D.; Hindmarsh, A. C. ACM Trans. Math. Soft. 1975, 1, 71.
(9)
where e is the elementary protonic charge, 0 is the permittivity of the vacuum, and D is the relative dielectric constant of the solvent water. Due to the differences in the Bjerrum length and the molecular models, ξ values in this paper are slightly different from those in the experimental reports. In excess salt cases, the results are presented as a function of X1/2, where X ) CP/CS. Although the unit of the concentration used in several experiments is given by the equivalent per 1 kg of water, the difference with our definition is ignored. Results The sodium ion (Na+) activities in several polyelectrolyte solutions without excess salt were measured by Nagasawa and Kagawa using the sodium amalgam electrode and the ion exchange membrane electrode. They pointed out the fact that the Na+ activity coefficient decreases with dilution of sample solution and it is entirely independent of the degree of polymerization.23,24 The counterion activity (23) Nagasawa, M.; Kagawa, I. J. Polym. Sci. 1957, 25, 61. (24) Rice, S. A.; Nagasawa, M. Polyelectrolyte Solutions: A Theoretical Introduction; Academic Press: London, 1961.
4126 Langmuir, Vol. 15, No. 12, 1999
Figure 6. Comparison with the data of Na-pectate-NaCl by Joshi and Kwak. Corrected activity coefficients in Na-pectateNaCl solution: γ-c, γ(c (closed symbols), and γ+c, from top to bottom. Small symbols, experiments; large symbols, MC-h; lines, Poisson-Boltzmann (a ) 3.0 Å) results. Solid lines and circles, CP ) 0.01 N; broken lines and triangles, CP ) 0.025 N. The DS value is given to be 1.00 (ξ ) 1.603) in the calculations.
Figure 7. Comparison of single ion activity coefficients of the MC-c and MC-h results with the Poisson-Boltzmann solutions in Na-CMC-NaCl. Solution conditions are the same as those in Figure 5. Closed symbols, MC-c; open symbols, MC-h; lines, Poisson-Boltzmann (a ) 2.5, 3.5, 4.5 Å, from bottom to top) results. Solid lines and circles, CP ) 0.01 N; broken lines and triangles: CP ) 0.025 N.
coefficients estimated from the Monte Carlo simulation are compared directly with the Na+ activity coefficient in the salt-free Na-CMC solutions. As shown in Figure 3, the Monte Carlo results agree sufficiently with the experiments not only on the absolute values but also on the concentration dependence. The apparent deviation of
Nishio and Minakata
the MC-h result from that of MC-c is not found in this salt-free case. The analytical Poisson-Boltzmann solution agrees with the Monte Carlo, especially MC-c, result most satisfactorily when the polyion radius, a, is chosen to be 3.5 Å. This polyion radius is slightly smaller than that of the molecular model in the Monte Carlo method. It seems that the agreement of the Monte Carlo and Poisson-Boltzmann results with the experiment is insufficient for the polyion having a flexible polymer skeleton.25 One of the available measurements of the single ion activities in the salt-containing CMC solutions was reported by Rinaudo and Milas. They obtained each single ion activity using the ion-selective electrodes for two CMC samples with different degrees of substitution (DS ) 1 and 2.5).26 The Monte Carlo and Poisson-Boltzmann calculations are attempted for their experiments under the same conditions of the NaCl solution. The dependence of the Poisson-Boltzmann solution on the polyion radius, a, is small. Here, the MC-h results agree with the PoissonBoltzmann solutions better than MC-c. The comparison of the MC-h and Poisson-Boltzmann results with their corrected activity data is shown for two samples in the presence of NaCl in Figure 4. The agreement of the mean activity is satisfactory for both samples with different degrees of substitution. However, the calculated single ion activity coefficients of the counterion and co-ion largely deviate from the experiments. The deviations of both ion activities cancel out in the mean activities. Joshi and Kwak presented a fine measurement for the Na-polyelectrolyte-NaCl system. They obtained mean and counterion activities using the ion-selective electrodes at two polyion concentrations for each sample.27 The MC-h and Poisson-Boltzmann calculations are shown with their data of Na-CMC-NaCl and Na-pectate-NaCl in Figures 5 and 6, respectively. The dependence on the polyion concentration is insignificant in both the experiment and the calculation. The calculated results of Na-CMC-NaCl are considerably close to the experiments for the coefficients of the mean activity and the counterion activity. However, the co-ion activity coefficients are overestimated in the simulations. Similar tendencies are observed in the case of Na-pectate-NaCl. The Poisson-Boltzmann solution at the polyion radius a ) 3.0 Å agrees with the Monte Carlo result well for the molecule model of the pectate. Under this relatively high CP condition, the surface effect in the cylindrical system appears more clearly in the low X range. Comparisons of single ion activity coefficients of the MC-c and MC-h results with the Poisson-Boltzmann solutions are represented for Na-CMC-NaCl in Figure 7. The deviation in the MC-c system is due to the surface effect on the ion distribution in the isolated cell system. The decreases in the ion concentrations are apparently observed near the cell surface of the cylindrical system at low X as shown in Figure 8a. At high X, the effect is not observed (Figure 8b). Discussion The present Monte Carlo study revealed that the mean activity coefficients as well as the counterion activity coefficients are in good agreement with several experimental results referred for the rodlike polyelectrolytes with relatively high linear charge density. This means (25) Nishio, T.; Minakata, A. Rep. Prog. Polymer Phys. Jpn. 1998, 41, 141. (26) Rinaudo, M.; Milas, M. Chem. Phys. Lett. 1976, 41, 456. (27) Joshi, Y. M.; Kwak, J. C. T. J. Phys. Chem. 1979, 83, 1978.
Ion Activities in Rodlike Polyelectrolyte Solutions
Figure 8. Comparison of ion distributions in the MC-c and MC-h systems with the Poisson-Boltzmann solutions for NaCMC-NaCl solution at CP ) 0.01 N. Other conditions are the same as those in Figure 5. (a) Local concentrations of counterion and co-ion as a function of the distance from the cell axis r at X1/2 ) 0.400. (b) The same at X1/2 ) 1.977. Small symbols, MC-c; large symbols, MC-h; lines, Poisson-Boltzmann (a ) 2.5, 3.5, and 4.5 Å from left to right) results. Crosses, circles, and upper lines, counterion molar concentration; plus signs, squares, and lower lines, co-ion molar concentration.
that the assumptions introduced in the canonical Monte Carlo calculation such as hard core polyion and small ions with definite sizes are fairly reasonable ones without losing fundamental reality in the polyelectrolyte solutions. Indeed, the potentiometric titration of CMC was well reproduced in a similar system.6 The results of the Monte Carlo and Poisson-Boltzmann calculations of the ion activity coefficients clarified that they depend significantly on the polyion charge density but slightly on the polyion radius as well as the polyion concentration under the conditions of the present study. Both calculations give almost the same results of the local ion concentrations at the cell boundary, if the surface effect can be eliminated in the Monte Carlo simulation and a suitable polyion radius is chosen in the Poisson-Boltzmann calculation. This implies that the difference in basic assumptions between both methods is not so significant in the present study. The Monte Carlo method gives a little different result of the concentration profile from that of the Poisson-Boltzmann calculation. This seems to be due mainly to the difference of the Monte Carlo system from the classical Poisson-Boltzmann approximation, such as the effect of the ion-ion correlation and of the finite ion size. To incorporate these effects, a modified Poisson-Boltzmann approach must be efficient.28 (28) Das, T.; Bratko, D.; Bhuiyan, L. B.; Outhwaite, C. W. J. Phys. Chem. 1995, 99, 410.
Langmuir, Vol. 15, No. 12, 1999 4127
Figure 9. Comparison of the Poisson-Boltzmann (a ) 3.5 Å) result with other theories. Corrected activity coefficients of CMC solution with mono-monovalent salt in the case of DS ) 1 (ξ ) 1.377) at CP ) 0.002 N. PB (solid line), Poisson-Boltzmann (a ) 3.5 Å) result; MLL (long-dashed line), Manning’s limiting law; NLL (short-dashed line), new limiting law by Iwasa et al.; GW (dash-dotted line), theory of Gue´ron and Weisbuch. Symbols at the end +, -, and ( indicate γ+c, γ-c, and γ(c, respectively.
The present results for the counterion activity are close to the new limiting law by Iwasa et al., although the coion activity coefficient is rather deviated.29 A little larger difference is found from the Manning’s limiting law.3 In the approximate formulas of the Poisson-Boltzmann equation by Gue´ron and Weisbuch, the co-ion activity coefficients are near to our numerical Poisson-Boltzmann results. Nevertheless, their counterion activity coefficient seems to be overestimated.8 Comparisons among these theoretical predictions at ξ ) 1.377 are summarized in Figure 9. As mentioned above, the results of the calculations well reproduce several experiments on Na-CMC and Napectate solutions for the mean activity coefficient and/or the counterion activity coefficients, although some unlike experimental data exist. In particular, the value of the co-ion activity is, however, overestimated in all cases. Similar observations are found in the modified PoissonBoltzmann calculations with the empirical correction in the case of the rodlike polyelectrolytes.30 In many studies, the co-ion activity coefficients were observed to be slightly below unity and were assumed to be almost unity in various theoretical works. This assumption may be adequate when CS is greater than CP. However, the co-ion activity coefficient larger than unity is obtained in our present framework when CS e CP. A similar situation provides an insufficient result for the simulation of the potentiometric titration of poly(glutamic acid).31 (29) Iwasa, K.; McQuarrie, D. A.; Kwak, J. C. T. J. Phys. Chem. 1978, 82, 1979. (30) Delville, A. Biophys. Chem. 1984, 19, 183.
4128 Langmuir, Vol. 15, No. 12, 1999
Concerning the co-ion activity, there are two points to be considered carefully. First, the experimental difficulty and ambiguity should be pointed out especially for the measurement of the single ion activity. Indeed, the corrected activity coefficients of the co-ion in the measurements by Joshi and Kwak except for CMC are almost unity. The second possibility is the insufficiency of the correction for the experimental data, expressed in eq 6. Iwasa and Kwak have proposed a more adequate method for correction based on the counterion condensation hypothesis.32 Nevertheless, the method is not fully sufficient to interpret the deviation of the co-ion activity.27 Reexamination for the correction of the co-ion activity surrounded by a large excess of counterions seems to be necessary for the experimental data.33 In conclusion, the present Monte Carlo study predicts fairly reasonable values for mean and counterion activity coefficients, but some reconsideration may be required (31) Nishio, T. Biophys. Chem. 1998, 71, 173. (32) Iwasa, K.; Kwak, J. C. T. J. Phys. Chem. 1976, 80, 215. (33) Alexandrowicz, Z. J. Polym. Sci. 1962, 56, 97.
Nishio and Minakata
for the co-ion activity coefficient either from the experimental accuracy or from the method of correction. If there is another possibility for the cause of the deviation of the co-ion activity, this may be related to the formation of the multichain domain at no or low ionic strength.34-36 When the attractive interaction between polyions suggested recently is clarified, the appropriateness of the cell model must be reconsidered.37 Acknowledgment. The authors are grateful to Professor J. C. T. Kwak and Dr. A. Rodenhiser for supplying the results of the measurements. The Monte Carlo calculations were executed on the vector processor in the information processing center of Hamamatsu University School of Medicine. LA981188R (34) Sedla´k, M.; Amis, E. J. J. Chem. Phys. 1992, 96, 817. (35) Sedla´k, M.; Amis, E. J. J. Chem. Phys. 1992, 96, 826. (36) Sedla´k, M. J. Chem. Phys. 1996, 105, 10123. (37) Ray, J.; Manning, G. S. Langmuir 1994, 10, 2450.
