J. Phys. Chem. 1992,96,1921-1932
the case whether continuum or molecular descriptions of the solvent were used. However, MD simulations of tryptophan using a united atom representation of hydrogen suggest the exact opposite, that x, interconversions are more frequent than x2 interconversions. The united atom representation of hydrogen was introduced into packages such as CHARMm and GROMOS because it significantly reduces the computational requirements of most problems pertaining to biomacromolecules. This model has been shown to present a satisfactory representation of certain properties such as structural minima. But empirical potentials are parametrized to reproduce properties near local/global minima; there is no guarantee that they will reflect the true nature of the potential surface far from local minima. Thus, there can be relatively good agreement between the explicit and the united atom descriptions of hydrogen, viz. the (x,, x2) geometries of the six rotamers of tryptophan, yet disagreement in regions separating these minima. The nature of the potential surface separating minima strongly influences trajectories of molecular dynamics simulations. Yet it is only with the advent of greater computer resources that some discrepancies, such as that illustrated here, may be exposed by sufficiently long simulations. Simulation times for proteins typically are of the order of tens of picoseconds, much shorter than the 1-ns simulations on the single tryptophan molecule described here. Even so, our simulations were not long enough to establish,
1921
with any statistical reliability, transition rates between rotamers on our potential energy surface. It therefore behooves users of M D simulations to carefully examine their results to establish consistency with experimental data whenever available. Otherwise, completely wrong interpretation and microscopic modeling of the experiments could ensue. For example, in this work, computational results using the CHARMm intramolecular potential (and full hydrogen representation) for tryptophan agree with the xi rotamer model used to interpret some N M R data.33 However, as found p r e v i o u ~ l y ,there ~ ~ is not even qualitative agreement between computational and experimental results concerning the predicted relative rotamer populations. Potential energy surfaces and M D simulations can offer experimentalists a rationale upon which to extend interpretations of their data. For example, it is conventional to interpret N M R coupling constants as arising from rotamers having one fixed geometry, as was done in ref 29. Figure 4 emphasizes that contributions to such experimentally measurable quantities actually come from a distribution of conformer geometries. Moreover, these distributions may not be narrow or symmetric (as is the case for perp-x2 g and anti-x2 t rotamers). Therefore, theoretical calculations can provide a framework upon which data analyses are developed that do incorporate a treatment of the librational motion that contributes to NMR relaxation times and fluorescence anisotropies.
Simulations on the Counterion and Solvent Dlstribution Functions around Two Simple Models of a Polyelectrolyte Heather L. Gordont and Saul Goldman* Department of Chemistry and Biochemistry, and the Guelph- Waterloo Center for Graduate Work in Chemistry, University of Guelph, Guelph, Ontario, Canada, Nl G 2 WI (Received: August 12, 1991)
We present our results from a series of long canonical ensemble Monte Carlo simulations on the counterion and solvent distribution functions in regions close to a high charge density negatively charged polyelectrolyte. The simulations are all for the "salt free" solution; i.e., a neutralizing number of counterions are present, but there is no additional supporting electrolyte. We considered two models for the solvent-SPC water and a dielectric continuum whose static dielectric constant equals that of SPC water-and two models for the polyion-a uniform line of charge along the central axis of a soft rigid cylinder and a helical necklace of charges spiraling around the surface of a soft cylinder. These models were used to carry out three sets of simulations: continuum solvent with axially charged polyion, discrete solvent with axially charged polyion, and discrete solvent with polyion containing the helical necklace of negative charges. Our results with the continuum solvent were consistent with those from a number of theoretical and simulation studies done in the past with similar models. This model produced a high narrow peak for the counterion concentration close to the surface of the polyion, which decayed rapidly with distance from the polyion. On the other hand, our simulations using (molecular) SPC water with a uniformly charged polyion p r o d u d a surprise. This time the counterions avoided regions near the polyion surface. This result was attributed to an unphysically high degree of solvent orientational polarization in this system. The highly polarized solvent particles successfully displaced all the counterions from regions near the polyion. This unphysical degree of solvent polarization turns out to be due to the use of a uniform line of charge to represent the charge distribution of the polyion. This was proven by our final set of simulations, where the charge distribution of the polyion was taken to be the spiraling helical necklace of discrete charges. In these runs the counterions were found to cluster near the polyion surface, and the extent of orientational polarization of the solvent was very much reduced from what was found using the axially charged polyion.
1. Introduction A polyelectrolyte is a macromolecule with regularly spaced ionizable groups. The purpose of this study is to learn about how a rigid, linear polyelectrolyte, in a water-like solvent, influences the distribution of solvent molecules and counterions near its surface. As with previous work on polyelectrolyte^,'-^^ the incentive to do this stems partly from the problem itself, which is 'Present address: Bioinformatics Section, Institute For Biological Sciences, Building M54,Montreal Road Campus, National Research Council of Canada, Ottawa, Ontario, Canada, K1A OR6.
intrinsically interesting, and partly from the DNA connection. Double-stranded DNA in water is a high charge-density polyelectrolyte, whose conformational, thermodynamic, and binding properties are influenced by the presence and distribution of dissolved electrolyte. This is a large field to which a number of groups have made contributions. To put our own interest in context we provide a brief outline of the history of previous work on this problem. The emphasis and direction taken by previous workers is a reflection both of their scientific taste, and the period during which the work was done. The earliest theoretical studies of which we
0022-3654/92/2096- 1921$03.00/0 0 1992 American Chemical Society
1922 The Journal of Physical Chemistry, Vol. 96, No. 4, 1992 are aware involved the application of the Poisson-Boltzmann theory to this problem.l,2 In these studies the polyion was taken to be a hard cylindrical rod of finite radius, with uniform surface charge density, and the solvent was a dielectric continuum. The counterions (small ions) were taken to be point charges that interacted with the instantaneous field of the polyion, but interacted with only a mean field with the other counterions. The mean field was obtained from the average charge density of the counterions. A concentric cylindrical boundary enclosed the polyion rod and small ions. These calculations resulted in a smooth continuous distribution of small ions about the polyion. Later, more refined Poisson-Boltzmann treatments were carried o ~ t . ~ ? In these studies, the small ions were given a finite size, and the small ionsmall ion interaction was approximated by using a fluctuation potential obtained from the self-atmosphere of t& small ions. These more refined treatments produced a more compact double layer; i.e., the counterion concentration at the polyion surface was enhanced and it then dropped off more rapidly with distance than was the case with the earlier, simpler, Poisson-Boltzmann treatments. The counterion condensation theory put forward by Manr~ing,~-lO seems to have been provoked by a comment communicated to Manning by Onsager.'l Onsager pointed out that, if the charge density of an infinite line of charge was high enough to make f > 1 (5 is a dimensionless measure of the polyion linear charge density), and if this line of charge was surrounded by point-charge counterions, the phase integral of the system would become singular, so that the system would be unstable. In order to eliminate this (physically impossible) instability, Manning postulated that sufficient numbers of (monovalent) counterions condense on the polyion so as to reduce the effective value of f (now comprising contributions from both polyion and condensed counterion charges) to unity, thereby removing the instability. Thus, counterion condensation theory results in a twestate model, in which a fraction ON (per polymer charge) of the counterions are 'bound" to the polyion and constitute a condensed phase. The remaining "free" counterions interact with the effective or net residual charge on the polyion. For polyions modeled as impenetrable cylinders containing a uniform axial line charge, Manning's theory is in the form of a (double) limiting law which predicts that, in the limit of infinite polyion length, the counterions will remain in the condensed phase, in the limit of infinite dilution. Despite this unphysical prediction, Manning's theory is correct for a polyion in the double limit of infinite polyion length followed by infinite dilution. This was proven independently by Ma~Gillivray,'~*'~ who obtained Manning's result (condensation persisting at infinite dilution) directly from the solution of the Poisson-Boltzmann equation. The model of an infinitely long cylinder of finite radius, at infinite dilution, is mathematically equivalent to a line of ~ h a r g e . ~Thus ~ ~ ~ ~ Manning's result, from the basis of a statistical mechanical approach, and MacGillivray's, from the basis of the PoissonBoltzmann equation, both supported the validity of Onsager's comment. The counterion condensation theory does not describe the distribution of the counterions in the condensed phase. The condensed fraction ON of counterions is solely a function o f f and the valence N of the counterions and is invariant with ionic strength. According to this theory, for cylindrical polyion models, (1) Fuoss, R. M.; Katchalsky, A.; Lifson, S. Proc. Nafl. Acad. Sci. U.S.A. 1951, 37, 579.
(2) Alfrey, Jr., T.; Berg, P. W.; Morawetz, H. J . Polym. Sci. 1951, 7 , 543. (3) Bratko, D.; Vlachy, V. Chem. Phys. Lett. 1982, 90, 434. (4) Outhwaite, C. W. J . Chem. SOC.,Foraday Trans. 2 1986, 82, 789. (5) Manning, G.S. J . Chem. Phys. 1969, 51, 924. (6) Manning, G . S. J . Chem. Phys. 1969, 51, 934. (7) Manning, G.S. J . Chem. Phys. 1969, 51, 3249. (8) Manning, G. S. Annu. Reo. Phys. Chem. 1972, 23, 117. (9) Manning, G. S. Q.Reo. Biophys. 1978, I f , 179. (10) Manning, G. S. Acc. Chem. Res. 1979, 12, 443. (1 1) Onsager, L. Footnote 13 of ref 5. (12) MacGillivray, A. D. J . Chem. Phys. 1972, 56, 80. (13) MacGillivray, A. D. J . Chem. Phys. 1972, 56, 83.
Gordon and Goldman the volume (per polyion charge) occupied by the condensed counterions is also salt concentration-invariant. However, Poisson-Boltzmann and other statistical mechanical treatments (see below) do provide descriptions of the counterion distribution function and find that the condensed volume does vary with the quantity of added salt. Recently the hypernetted chain (HNC) integral equation has been applied to a cell model of a dilute polyelectrolytesolution.'e19 The principal advantage of this approach is that it provides a better treatment of the ion-ion correlations than is possible with the Poisson-Boltzmann equation. The model for the polyion is again ~an infinitely long, rigid, uniformly charged hard or soft rod, surrounded by small ions at a finite concentration in a dielectric continuum. The small ions are hard or soft spheres, each with an imbedded point charge. The small ionsmall ion correlations were treated both by the (simpler) mean spherical approximation (MSA)1618 and by the full HNC appr~ximation.'~.'~,~~ The small ion-polyion correlations were always treated with the full H N C approximation. These two treatments agree qualitatively with each other, and with the more sophisticated of the PoissonBoltzmann treatments,3q4in that they all predict that the double layer around the surface of the polyion is more compact than the one predicted by the simplest Poisson-Boltzmann-based calculations.1*2The latter calculations predict too low a concentration of counterions (by 15%) adjacent to the polyion surface. The HNC results predict a continuous distribution of counterions about the polyion, as opposed to a condensed phase or two-state distribution. These results also concur with those from the Poisson-Boltzmann equation that there is no discontinuity in the counterion distribution function around f = 1, and that the fraction of counterions near the surface of the polyion increases continuously with the concentration of added salt. Previous simulation studies on polyelectrolytes are of two types. In most of these studies the solvent was modeled as a dielectric c o n t i n ~ u m , ~ but J ~ -in~ a~ a discrete molecular model of water was used. The continuum model simulations were performed mostly to check the accuracy of the previously mentioned theories. Most used an infinitely long, rigid cylinder with a uniform linear charge along the axis as the model for the polyion, which was surrounded by a cylindrical annulus that contained a finite concentration of small ions. The results of these simulations support the results on the small ion distribution obtained from either the modified Poisson-Boltzmann calculations or the H N C calculations. Two grand canonical Monte Carlo (GCMC) simulations on the continuum solvent model have been reported?2.24from which small ion activity coefficients and Donnan exclusion coefficients can be obtained. The results of these GCMC simulations supported those from the Poisson-Boltzmann and the HNC theories. ~ 'There ~ have also been several simulations, in the primitive solvent (dielectric continuum) regime, designed to test the validity of a
-
(14) Baquet, R. J.; Rossky, P. J. J . Phys. Chem. 1984, 88, 2660. (15) Baquet, R. J.; Rossky, P. J . J . Phys. Chem. 1988, 92, 3604. (16) Lozada-Cassou, M. J . Phys. Chem. 1983, 87, 3729. (17) Gonzales-Tovar, E.; Lozada-Cassou, M.; Henderson, D. J . Chem. Phys. 1985, 83, 361. (18) Vlachy, V.; McQuarrie, D. A . J . Chem. Phys. 1985, 83, 1927. (19) Murthy, C. S.;Baquet, R. J.; Rossky, P. J. J . Phys. Chem. 1985,89, 701.
(20) Mills, P.; Anderson, C. F.; Record, Jr., M. T. J . Phys. Chem. 1985, 89, 3984. (21) Mills, P.; Paulson, M . P.; Anderson, C. F.; Record, Jr., M. T. Chem. Phys. Lett. 1986, 129, 155. (22) Mills, P.; Anderson, C. F.; Record, Jr., M. T. J . Phys. Chem. 1986, 90, 6541. (23) Paulsen, M. D.; Richey, B.; Anderson, C. F.; Record, Jr., M. T. Chem. Phys. Lett. 1987, 139, 448. Erratum. Chem. Phys. Lett. 1988, 143, 115. (24) Vlachy, V.; Haymet, A. D. J . J . Chem. Phys. 1986, 84, 5874. (25) Nilsson, L. G.;Nordenskiold, L.; Stilbs, P. J . Phys. Chem. 1987, 91, 6210. (26) Seibel, G.L.; Singh, U . C.; Kollman, P. A. Proc. Narl. Acad. Sci., U.S.A. 1985, 82, 6531. (27) Forester, T. R.; McDonald, I. R. Mol. Phys. 1991, 72, 643. (28) van Gunsteren, W . F.; Berendsen, H. J. C.; Geurtsen, R. G.; Zwinderman, H. R. J . Ann. N . Y . Acad. Sci. 1986, 482, 287.
Counterion and Solvent Distribution Functions
region containing water and/or small ions
I
\L -Z
Figure 1. Schematic diagram of simulation cell.
uniform line of charge as a representation of the polyion’s Coulombic field.21925929J0Of these, only Mills et ala?*whose helical model was taken to represent the B form of DNA, and Nilsson et al.,z5 who modeled poly(styrenesulfonate), correctly incorporated long-range electrostatic interactions. Both groups found that the concentration of univalent counterions, close to the polyion surface, increased moderately over what was found for the line charge model. The relative enhancement was greater, the greater the charge density on the polyion. There is, a t present, substantial disagreement concerning the nature of the counterion distribution function, among molecular dynamics simulations on models of DNA surrounded by a discrete water-like solvent.2628 The models of the polyion used in these simulations are of the site-site type and are intended to model particular forms of DNA in water, as opposed to polyions in general. Specifically, Seibel et a1.26report a much stronger degree of attachment of a Na+ ion to a phosphate group oxygen than is reported by Forester and M~Donald.~’At the other extreme, van Gunsteren et a1.28report that Na+ ions always remain solvent-separated from phosphate groups. While these three studies have the Na+ ions hovering around the polyion, the mean distance of a Na+ ion from a phosphate oxygen (for example) ranges from tight contact to solvent-separated, depending on the study and the details of the model. Our interest is to understand, in a general way, how a polyion influences the distribution of the surrounding solvent and counterions. To do this, we set up a series of simulations in which we systematically vary first the model for the solvent and then the model for the polyion. In one set of simulations we fix the model for the polyion as an infiiite, rigid, soft cylinder with an imbedded axial line of charge, and use this with both a discrete solvent (rigid SPC water3’) and a dielectric continuum whose static dielectric constant is that of rigid SPC water (t N 6632-34). This series of calculations is done to determine how a discrete solvent and dielectric continuum differ in terms of their influence on the counterion distribution function close the surface of a polyion. A second set of simulations is designed to determine how a different polyion field that arises from a different charge distribution within the polyion affects the counterion and solvent distribution functions. In these calculations we use rigid SPC water as the solvent, and a helical arrangement of charges on a cylindrical surface, as the model for the polyion. A subsequent comparison of the results of this calculation with the one in which an axial line of charge was used for the polyion charge distribution gives the influence of the polyion charge distribution on the solute and solvent distribution functions. (29) LeBret, M.; Zimm, B. H . Biopolymers 1984, 23, 271. (30) Conrad, J.; Troll, M.; Zimm, B. H. Biopolymers 1988, 27, 171 1 . (31) Berendsen, H. G. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. Interaction models for water in relation to protein hydration. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, 198 1 . (32) Watanabe, K.; Ferrario, M.; Klein, M. L. J . Phys. Chem. 1988, 92, 819. (33) Gordon, H.; Goldman, S.Mol. Simuln. 1989, 2, 177. (34) Alper, H. E.; Levy, R. M. J . Chem. Phys. 1989, 91, 1242.
The Journal of Physical Chemistry, Vol. 96, No. 4, 1992 1923 2. Simulation Cell A schematic representation of the simulation cell is shown in Figure 1. The simulation cell is a cylinder whose outer wall is concentric with the cylindrical polyion; the central axis of the polyion is also the axis of the simulation cell. The small ions (counterions) and mobile SPC particles were confined to a cylindrical annulus around the polyion and could not escape through the outer cylindrical wall of the cell. We used two types of cylindrical outer walls, depending on whether the solvent was a dielectric continuum or discrete SPC water. For the continuum solvent, the outer boundary was a smooth impenetrable wall. For the discrete solvent we used a ‘glassy wall” for the outer boundary. The glassy wall was a cylindrical annulus, concentric with the polyion and its surrounding solution, composed of immobile SPC particles, 3 solvent diameters (3uspc) thick. These SPC particles were mathematically carved out from a large (2420 particles) homogeneous cube of SPC water of density 1 g/cm3 that had been extensively equilibrated a t 300 K. The particles selected to comprise this outer annular wall were those whose centers were situated within the coordinates of the desired annulus. This wall had an inner radius of about 25 A, thickness 3uspc, and height 25.5 A. There were 1465 particles in our glassy wall. The positions and orientations of these particles were permanently fixed. The mobile SPC molecules and counterions surrounding the polyion were thus contained in the radial direction of the simulation cell by a cylindrical shell of SPC water, frozen in a configuration typical of bulk SPC water. Hence the term ‘glassy wall”. We note in passing that a glassy wall in the form of a spherical shell has been used previously in other Periodic boundary conditions were used in the z-direction; the z-axis pointed along the polyion axis. The polyion was immobile and trial steps in the Monte Carlo simulations consisted of selecting and moving counterions and, when the solvent consisted of discrete particles, mobile SPC particles. The z-coordinate of the selected particle was always zero, since the effect of using periodic boundary conditions in the axial direction was to reset the vertical center of the simulation cell to zero for every state sampled. Our results are for a simulation cell whose axial length (Lin Figure 1) was 25.5 A and whose radial distance between the axis of the polyion and the constraining wall (Roin Figure 1) was exactly 25 A for the smooth wall, and approximately 25 A for the glassy wall. Since the polyion had a radius of 10 A (below), the solution around the lyion was contained in a cylindrical annulus of height 25.5 and thickness 15 A.
sf“
-
-
3. Models (a) Polyion. Two models for the polyion were used. Both were infinitely long, had finite radii (- 10 A), and had soft, cylindrical Lennard-Jones-like cores. These cores revented solvent or counterions from penetrating beyond 2-3 into the body of the pol yion. In one model, the polyion’s negative charge was distributed uniformly along its central axis. The charge density was taken to be -e/6, where e is the electronic charge and 6 is the length of a segment of the axis that carries exactly one unit of charge. To mimic, crudely, the charge density of DNA, 6 was taken to be 1.7 In the second model, the total charge on the cylinder segment inside the simulation cell was contained in 15 fixed soft-core ions, which were distributed in the form of a helical necklace that spiraled around the outside cylindrical surface of the polyion. Each ion had a charge of -e and the center of each ion was 10 A from the polyion axis. In this second model the polyion segments that stretched from the vertical extremities of the simulation cell to +m or -m were the same as in the first polyion; Le., they were smooth cylinders with a uniform line of charge along the central axis. Thus our two polyion models differed from one another in two respects, one major and the other minor. The major difference was in how the polyion charge was distributed. The minor dif-
R
A.19320
(35) Sheykhet, L.; Simkin, B. J . Mol. Liq. 1988, 37, 153.
1924 The Journal of Physical Chemistry, Vol. 96, No. 4, 1992
Gordon and Goldman
Therefore, we used a glassy wall for our molecular solvent calculations, both to contain the system radially and to reduce these eLJ, erg ULJ, A structural artifacts. It turned out that the glassy wall very much 108 X 3.166 particle-particle* reduced these artifacts, relative to what would have been obtained 18.526 X 11.583 particle-polyionc with a smooth wall. This improvement has been documented in an earlier p~blication.~’We did use a smooth constraining wall “‘Particle” means SPC water or counterion. bThese are the Lenfor our dielectric continuum calculations, because wall-induced nard-Jones core parameters of the SPC potential. They were used as the Lennard-Jones part of the counterion-counterion, SPC-SPC, artifacts (in a continuum solvent), of course, do not arise. counterion-SPC, discrete polyion charge-counterion,and discrete poThe glassy wall is an approximate representation of the discrete lyion charge-SPC potential energy expressions. These values were solvent extending to infinity in the radial direction. As with other selected so that the ion-polyion potential minimum occurred at apapproximations, such as smooth hard walls or periodic boundary proximately 10 A. This was done to facilitate comparisons with other conditions, this approximation is bound to have some influence work (see text). on the results. There are two reasons for our view that the influence of the glassy wall approximation on the counterion disference was that the uniformly charged polyion’s surface was tribution function is relatively unimportant. First, we carried out everywhere smooth, whereas that of the polyion with the helical a series of calculations to determine the role of the size of the charge distribution had a bumpy surface inside the simulation cell. simulation cell on the counterion distribution function. In these The bumps were the protruding parts of the Lennard-Jones spheres calculations, we used a continuum solvent, the salt-free solution, that were centered on the surface of the polyion (10A from the a smooth hard wall to contain the system in the radial direction, central axis) which contained the embedded negative point charges. and periodic boundary conditions in the axial direction. We found (b) Solvent. Two models for the solvent were used. One was that variations of the size of the radius (Roin Figure 1)had very a dielectric continuum whose static dielectric constant equalled small effects on the counterion distribution function. Specifically, that of rigid SPC water (66 at 25 OC and 1 g m / ~ m ~ ) . ~This *-~~ a reduction of Ro from 100 to 25 A, using e = 78 and L = 25.5 model was used only with the polyion whose charge distribution A, produced no discernible change in the shape of the distribution was a uniform line of charge along the central axis. function, no discernible change in the peak position, and a change The other model was rigid SPC water?’ which was used with of only 15% in the height of the counterion concentration peak both polyion models. (the smaller cell produced the higher peak). These results, shown (c) Countenom. All our calculations involved 15 positively in Figure 7.1 of ref 45,were taken to indicate that, provided Ro charged counterions, each with a charge +e. Since the charge was not less than -25 A (for L = 25.5 A), the counterion dison the polyion within the central cell was -15e, the use of 15 tribution function was not highly sensitive to the approximation univalent counterions ensured electroneutrality within the central used for the radial boundary. The second source of support comes cell. from the form we will adopt to present our results. Specifically, Each counterion consisted of a point charge (+e) embedded we will compare the counterion distribution obtained for two a t the center of a Lennard-Jones sphere. The values of the polyion models, for which the same glassy wall is used to contain Lennard-Jones E and u parameters were taken to be those of rigid the system radially. Thus, any artifacts due to the glassy wall SPC water (Table I). This was done to avoid specific packing will at least partially cancel in the comparison. or solvation effects, which might otherwise have obscured the (c) Simulations Using a Dielectric Continuum for the Solvent. primary effects we sought to elucidate. These calculations were canonical ensemble (fKed N, V, 7‘)Monte Carlo simulations, in which a variant of the Metropolis Monte 4. Details of the Methodology and Simulations Carlo (MMC) algorithm42 was used. The variant, which we In this section we complete the description of our models and developed and described in detail p r e v i o ~ s l yis , ~called ~ densityapproximations and provide the technical details that underlie our scaled Monte Carlo (DSMC). DSMC differs from MMC in that simulations. A number of these details have been described MMC uses a sampling cube of fixed size, while DSMC involves p r e v i ~ u s l yso~ that ~ , ~the ~ description below is brief. The detailed scaling the size of the sampling cube to the density of the region potential energy expressions we used are of interest to the specialist, being sampled. As with MMC, DSMC obeys detailed balance but not, perhaps, to the general reader. Therefore, these exand so ultimately converges to the correct limiting distribution. pressions are given in the Appendix to this article. However, DSMC converges much more rapidly than MMC for (a) The Model for the Helical Polyion. As seen from section simulations where the particle density varies significantly within 3a, our charge distribution was somewhat disjointed for this model. the simulation cell.38 In our application we used small sampling The charge distribution was a set of discrete charges on a helix cubes near the polyelectrolyte where the counterion density was inside the simulation cell, but was a continuous axial charge high, and large sampling cubes, in regions remote from the podistribution above and below the cell. We did this to simplify the lyelectrolyte surface, where the counterion density was low. counterion-external polyion and SPC-external polyion potential Specifically, the maximum move size in any Cartesian direction energy expressions. The approximation was justified by prelim(Ao) was given by inary trial calculations, from which we found that this simplification had only a very small effect on the total polyion field at 4 = (1/5)rip (1) the vertical center of the simulation cell, where it is relevant. where ri4is the distance of the counterion from the polyion axis. (b) The Glassy Wall Bodary. Our system was inhomogeneous This choice produced average acceptance rates in the range 0.3 in that it contained only one polyion at the center of the simulation to 0.7 over the entire radial extent of the simulation cell. The cell. Consequently, periodic boundary conditions are inappropriate simulation cell was described in section 2, the static dielectric in the radial direction, and some sort of cylindrical constraining constant e in eqs A5,A13,and A14 was 66,and N , the number wall had to be used. For our calculations with the molecular of counterions in the cell, was taken as 15 to ensure electroneusolvent, a smooth constraining wall would have been ill-advised, trality. No additional ‘supporting electrolyte” was present, either because smooth surfaces are known to perturb both the radial and here or in the runs described below. angular distributions of water-like liquids near them.39-4’ We used a semirandom initial counterion distribution (i.e., random except for disallowed counterion-polyion overlaps), and (36) Torrie, G. M.; Valleau, J. P. J. Chem. Phys. 1980, 73, 5807. equilibrated the system for 1.5 X lo6 trial steps, prior to collecting (37) Gordon, H.; Goldman, S. Mol. Simuln.1989, 3, 315. data. The averaged internal charge density was used to generate (38) Gordon, H.; Goldman, S. Mol. Simuln.1989, 3, 213. (39) Lee, C. Y.; McCammmon, J. A.; Rossky, P. J. J . Chem. Phys. 1984, for use in eq A13. The averaged external the function pext(r’) TABLE I: Lennard-Jones Parametersa
80, 4448.
(40) Valleau, J. P.; Gardner, A. A. J . Chem. Phys. 1987, 86, 4162. (41) Torrie, G. M.; Kusalik, P. G.; Patey, G . N. J . Chem. Phys. 1988,88, 7826.
(42) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J . Chem. Phys. 1953, 21, 1087.
Counterion and Solvent Distribution Functions
The Journal of Physical Chemistry, Vol. 96, No. 4, 1992 1925
charge density of the equilibration run was used to construct the Run A over-screened starting pext(r3for the subsequent collection run. During both equilibration and collection runs put(‘’) was periodically updated 3.0 by setting it equal to the cumulatively averaged pi&’). The integration required in eq A13, at each step, was done numerically by two-dimensional compound Gauss quadrature. We used 10 Gauss points for the angular integration and 40 Gauss E U points for the radial integration. In all, 9 X lo6 trial steps were 5 00 run in the data collection phase of the simulation. Run B over-screened F 4 0 (a) Smulatim with SPC Water and a Polyim with a Uniform e Axial Charge Dedty. The simulation cell contained 1454 mobile g 30SPC particles, and 15 univalent positive counterions. The glassy +wall, which contained the system in the radial direction, contained 8 200 1465 immobilized SPC particles. The usual periodic boundary 1 0 conditions were applied in the axial direction. 0 The three Monte Carlo simulations described in this section c ki 0 0 : Run C under-screened were for a canonical ensemble and involved singleparticle moves with force and torque biasing,@ with the biasing partially turned on (A = The torque biasing was done by the method given in the appendix to ref 44. The forces and torques were calculated analytically from the relevant expressions for the potential (Appendix). Expressions for the forces and torques have been written out e l s e ~ h e r e . The ~ ~ magnitudes of the maximum translational and rotational step sizes were adjusted so that the average ac00 ‘h2 5 30 5 10 15 20 ceptance rates anywhere in the cell were in the range 0.3 to 0.7, r,p(Q both for the SPC and the counterion moves. Verlet neighbor lists& with a list sphere radius of 3 . 5 were ~ ~ ~ Figure 2. Counterion concentration as a function of distance ripfrom the polyion axis. (-) discrete SPC water (---) continuum solvent with e used to efficiently locate pairs of SPC particles located within the = 66. Thex are the unsmoothed results, obtained by averaging over an cutoff distance, r,, from each other. The value of r, was taken annular thickness of 0.05 A. as 2 . 9 4 (~u s~p ~~= 3.166 A). The integrations in eqs A13 and A22 were done by two-diand 1.75 X lo7 trial steps for runs A, B, and C, respectively. In mensional Gauss quadrature, using 32 Gauss points for the angular each case data collection was started after the functions we integration, and 5 sets of 32 Gauss points for the radial integration. moditored (section 5 ) stopped drifting. Additional details on this The runs were at 298.15 K and a t a solution density of approxseries of runs are given e l s e ~ h e r e . ~ ~ imately 1 g/cm3. This gross density was achieved by placing 1469 (e) Simulations with SPC Water and a Polyion with a Helical SPC particles semirandomly in the cylindrical annulus that Charge Distribution. The simulations here were very similar to comprised the cell. (It is impossible to be exact about the density those described above, except that here the potential energy exbecause our polyion and glassy wall were soft. This softness makes pressions for the internal interactions were t h w given in subsection the volume of the cell an imprecisely known quantity.) Then, 15 c of the Appendix rather than in subsection b. We did two runs, of the SPC particles were randomly picked and converted to both using the overscreened approximation for the external incounterions. teractions, which were dealt with exactly as in section 4d. The For clarity we label the three runs A, B, and C. Runs A and two runs differed in their starting configurations, which were taken B both were done using the “overscreened approximation” (Apfrom one configuration of equilibrated SPC water around a smooth pendix, subsection b), but differed from one another in their initial uncharged polyion, contained in the radial direction by our glassy configurations. Runs B and C were started from the same initial wall. For each run, a different set of 15 SPC particles, picked configuration, but differed in that run C was done using the at random, were converted to positive counterions, and 15 fixed “underscreened approximation” (Appendix, subsection b). The negative ions, which comprised our helix of charges, were disstarting configuration for run A was that of 1469 SPC particles tributed in accordance with eq A23 on the polyion surface. The preequilibrated at 300 K around the uncharged polyion, using the creation of the Lennard-Jones spheres around these negative glassy wall boundary. Then 15 of these particles were picked at charges resulted in overlaps with a few nearby mobile SPC random and converted to counterions. The starting configuration particles. These overlaps were slow to relax and so were eliminated for runs B and C was that of 1469 SPC particles equilibrated at prior to continuing with the equilibration by discarding the 300 K around the charged polyion, again using the glassy wall overlapping mobile SPC particles. Thus we proceeded with 1450 boundary, followed by randomly picking 15 particles for conversion SPC particles in one run and 1451 in the other (vs 1454 in the to counterions. Thus the solvent was started from a much more runs in section 4d). This difference (-0.3%) in gross solvent strongly polarized configuration in runs B and C than in run A. density relative to that in subsection d was considered unimportant.. Finally, we tried to speed up the convergence rate of the runs The runs, which we label D and E, were each equilibrated for by picking counterions for a trial move ten times more frequently lo7 trial moves and then data was collected for 7.5 X lo6 trial than SPC’s. We did this because the counterions made up only moves in each run. The acceptance rates, as before, were in the about 1% of the mobile particles around the polyion. Thus, 10% range 0.3 to 0.7 both for counterions and SPC particles. As in of the trial moves were made with counterions, and 90% with SPC section 4d, 10% of the trial steps were made with counterions. particles. Run A required a very long equilibration phase (36 X 106 trial 5. Results and Discussion steps) relative to runs B or C (13 X 106 trial steps each). Following Our results are given graphically in Figures 2-13. Figure 2 this equilibration phase, data were collected over lo7, 2 X lo7, shows the counterion distributions we obtained from our runs for a smooth polyion with a uniform axial charge distribution, for both a continuum and a molecular model of the solvent. Clearly, (43) Rao, M.; Pangali, C.; Berne, B. J. Mol. Phys. 1979 37, 1773. we got a very different result for the counterion distribution, (44) Goldman, S. J. Chem. Phys. 1983, 79, 3938. (45) Gordon, H. L. Monte Carlo Simulations of a Model Polyelectrolyte depending on the model used for the solvent. Using a Discrete Molecular Solvent. Ph.D. Thesis, University of Guelph, Our results using the dielectric continuum are qualitatively in 1990. accord with those of previous w o r k e r ~ . ~ . ~ , ’ ~ , ~Specifically, ~,l*-*~ (46) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; the counterions are found to cluster close to the surface of the Clarendon Press: Oxford, U.K., 1987; Section 5.3.1. L
v
r
’
s
Gordon and Goldman
1926 The Journal of Physical Chemistry, Vol. 96, No. 4, 1992 3.0 Run A:
over-screened
polyion surface
Ro
Run E: over-screened
Figure 5. Examples of (hypothetical) orientations of SPC water mole-
n
Run
C: under-screened
Figure 3. Singlet distribution function of SPC water as a function of distance from the plyion axis. If p is the bulk density then pdl)(rwp) is the &"pic average density of SPC particles at a distance rwpfrom the polyion axis. These curves are for averages over an annular thickness of 0.05 A and are not smoothed.
cules relative to n, the vector normal to the polyion surface. e,, is the angle between n and the dipole vector 1;Bo, is the angle between n and each of the vectors pointing along the OH bonds.
counterions, which are forced away from the polyion and toward the glassy wall. This explanation is supported by the structural and energetic results displayed in Figures 3-9. Figures 3-7 show the radial and angular distribution functions for SPC water in the simulation cell. Figure 3, the singlet distribution function of SPC water, shows that the uniformly charged polyion very tightly holds SPC water near its surface. This is best appreciated by comparing the results in Figure 3 with thase shown in Figure 4, where the radial structure of SPC water around an uncharged polyion (i.e., a soft cylinder) is shown. The high narrow fmt peak seen for the charged polyion vanishes when the polyion charge is turned off. The angular distribution functions that characterize the mean orientation of the SPC particles are shown close to, and far away from, the polyion in Figures 6 and 7, respectively. We use two orientational distributions functions, P,(B) and PoH(e).39'41These are normalized through 1 = l K P u ( e )sin 8 dB
1.0
(2)
3
L v
2 = J*PoH(e) sin
A r-
v
0
0.5
polyion charge off
~
e dB
and are defined by (3) i
Figure 4. Singlet distribution function of SPC water as a function of distance from the polyion axis for an uncharged polyion. The results are averages over an annular thickness of 0.1 A. Additional details on these simulations are given in ref 37.
polyion, creating a peak in the counterion concentration adjacent to the polyion surface. The concentration then drops off rapidly with distance from the polyion. The picture is of a concentrated compact double layer around the polyion that decays rapidly with distance. In contrast, the counterions were found to avoid regions near the polyion when a molecular solvent was used. This result was quite unexpected both on physical grounds, and because nothing like it has ever been reported. We are, however, certain that this result is correct for the model that produced it. It was obtained from long collection runs that followed long equilibration periods (section 4d) and it was obtained for both our underscreened and overscreened approximations and from very different starting configurations. Despite the fact that it is counterintuitive, this result is not difficult to understand. It is produced by an unrealistically high degree of solvent polarization, which in turn results from a physically unrealistic model for the polyion, namely, one that is uniformly charged. The highly polarized SPC particles monopolize the available space near the polyion, at the expense of the
where N(Oi) is the number of p or OH vectors making an angle of Bi & A4/2 with the surface normal. Here, A8 = IOo for all values of i. The distribution functions P J e ) and POH(e)should have constant values of 0.5 and 1.0, respectively, for a system displaying no angular anisotropy, such as a homogeneous liquid. The meaning of the angles and vectors is illustrated in Figure 5 . It is clear from Figure 6 that close to the polyion surface P,(B) and POH(B)have large peaks at 180° and 125O, respectively, which means that close to the polyion the dipole vector of SPC water tends to point directly toward the polyion axis (see upper left diagram in Figure 5). Figure 7 shows that sisnificant remnants of this orientational structure persist throughout the bulk of the solution. This is direct structural evidence that the uniform line charge model for a polyion very strongly polarizes SPC water. We have decomposed the total potential energy of an SPC particle into its various parts, and this breakdown is shown, for a part of one of our runs, in Figure 8. From this graph it is seen that the very low reduced potential energy of an SPC particle (plot 5) comes largely from the contribution of the SPC-polyion interaction (plot 4). Also, plot 3 of this figure shows that the SPC-SPC interaction weakens as the polyion is approached. This is due to the polyion polarizing the SPC particles, so that the SPC-SPC Coulombic interaction is significantly weakened relative to its mean value in the bulk liquid. An alternative and worthwhile way of discussing the foregoing result is in terms of the electrostatic potential throughout the simulation cell. This was obtained by storing a series of sto-
-
-
The Journal of Physical Chemistry, Vol. 96, No. 4, 1992 1927
Counterion and Solvent Distribution Functions
e e e e e F i p e 6. Angular distribution functions that characterizethe average orientations of SPC water molecules near the smooth, cylindrical, axially charged polyion surface. The functions PJ9) and POH(9) are defined by eqs 2 and 3 A9 = 10’. Legend: ( 0 )results using ‘underscreened” approximation, run C; (v)results using ‘overscreened” approximation, run B. In a homogeneous fluid, P,,(9) = 0.5 and POH(6) = 1.0 at all 9.
e e e e 0 Figure 7. Angular distribution functions that characterize the average orientations of SPC water molecules near the glassy wall boundary. Other details as in caption to Figure 6.
TABLE II: Pair Interactions Included in Our Three Sets of Calculations model polyion with uniform axial charge density and dielectric continuum solvent; equations in Appendix, section a polyion with uniform axial charge density and SPC water inside simulation cell; equations in Appendix, section b polyion with discrete helical charge distribution and SPC water inside simulation cell; equations in Appendix, section c
internal interactions’ mun terion-counterion coun terion-polyion
external interactionsb counterion-counterion charge density counterion-polyion
counterion-counterion counterion-polyion counterion-SPC SPC-pol yion SPC-SPC counterion-counterion counterion-polyion counterion-SPC SPC-polyion SPC-SPC
counterion-counterion charge density counterion-polyion SPC-counterion charge density SPC-polyion counterion-counterion charge density counterion-ply ion
SPC-counterion charge density SPC-polyion
“These entries, in the form a-b, mean that both parts of the pair a-b are inside the simulation cell. bThese entries, in the form a-b, mean that a is the part of the pair inside the cell, and b is the part outside (above and below) the simulation cell. chastically independent configurations 0 ( lo6 trial steps apart) during the collection phase of our runs, and then, for these configurations, working out the electrostatic potential a t the center of SPC particles at various locations in the cell. (This was done by switching locations of a counterion with all the SPC particles in the cell and working out the electrostatic part of the potential energy of interaction of this test ion with its environment. We then average this quantity at each of a series of radial distances from the polyion axis). As seen from the upper graph in Figure 9, the average electrostatic potential is very large and positive near
the polyion surface, falls with distance, and is negative near the glassy wall. This is consistent with the fact that, in this model, the counterions avoided the polyion in favor of the glassy wall. For comparison, we also show in the upper graph of Figure 9 the electrostatic potential for the continuum solvent model. It, of course, is negative near the polyion and decays toward zero with distance from the polyion, and its magnitude spans a much smaller range than was the case with the discrete solvent. In the lower graph of Figure 9 we show the breakdown of the average electrostatic potential in terms of the contributions by the polyion,
1928 The Journal of Physical Chemistry, Vol. 96, No. 4, 1992
Gordon and Goldman 10.0
R u n 0 : over-screened
[CI] 1
R u n E: over-screened
5
25
20
10
Figure 8. Components of average potential energy in reduced units of an SPC molecule as a function of radial distance from the polyion. Calculated for the 'under-screened" approximation, over 0.15 X 106 trial steps. Legend (1) SPC-extemal counterion charge density contribution; (2) SPC-counterioncontribution;(3) SPC-SPC contribution; (4) SPCtotal polyion contribution; ( 5 ) total SPC energy. 200
I
I
7
1
I
,
I
,
7
8
a
,
,
'
,
1
8.0
4.0 -
0.0
5
'*.
175 (De' 150 . ' e x i8 125 Y 75 . 50 . 25 .
30
polyion axis, using SPC water as the solvent, and a helical charge distribution on the surface of a soft cylinder, as the model for the polyion. Runs D and E had different initial configurations (see text).
'Y,
.*--
1.5
I
I
R u n D:
. "
25
F i e 10. Counterion concentration as a function of distance from the
0 .
-25 -50
20
15
10
rip6)
\
100.
-
2.0
"
'
'
"
.
"
counterion
.
.
'
"
over-screened
A
'
i
0'5
I i R u n E: over-screened
polyion - 1 0 0 0 ~ ' ' " ' ~ " ' ' ' 9 10 1 1 12 13 14 15 16 17 18 19 20 21 22 23 24 25
R(A)
Figure 9. The electrostatic potential as a function of distance from the
polyion axis, at the vertical center of the simulation cell, for the polyion with the linear charge density. Upper figure: ( 0 )Total electrostatic potential using the molecular solvent; (-) total electrostatic potential for a continuum solvent with a static dielectric constant of 66. Lower figure: breakdown of contributionsto the total electrostatic potential using the molecular solvent. ( 0 )Total electrostatic potential. the counterions, and the solvent. Note the very large positive contribution from the solvent near the polyion, which is a manifestation of extensive solvent polarization. This contribution, together with that of the counterions, is sufficient to offset the negative potential of polyion and thereby produces a net positive electrostatic potential near the polyion. We guessed that the extensive polarization of the solvent, which ultimately led to the counterions avoiding the polyion surface, was a consequence of our use of a uniform axial charge density as the representation of the charged part of our polyion. We used this model because it, or its equivalent, a uniformly charged cylinder, had been extensively used by other workers in the past.'-24 However, these previous studies all were for the primitive electrolyte model in which the solvent is replaced by a dielectric
0'5
I i
0.0
5
10
15
20
25
30
rwp(A) Figure 11. Singlet distribution function of SPC water around polyion
with helical charge distribution. continuum. Of course, for a continuum solvent polarization cannot arise, so that even a crude representation of the polyion's charge can still produce a qualitatively correct result for the counterion distribution function. To demonstrate that this was in fact the case, we carried out two additional simulations, both with SPC water as the solvent, but this time with a semirealistic model for the charge distribution on the polyion, namely the helical necklace of charges described in subsection c of the Appendix and sections 3a, 4a, and 4e. The results of these simulations are shown in Figures 10- 13. It is seen from Figure 10 that, for the helical necklace model of the polyion, with a molecular solvent, the counterions cluster close to the surface of the polyion and form a concentrated double layer near it. The small peak at ~ 2 A4is due to one of our 15
The Journal of Physical Chemistry, Vol. 96, No. 4, 1992 1929
Counterion and Solvent Distribution Functions
, 9.0(rwp(l O.OA
8.0(rwp(9.OA
-30
'
8
1O.O(rwp(l1 .OA
I
10
12
14
16
1000 I
18
20
22
24
26
1
counterion
Figure 13. Average electrostatic potential ih the simulation cell for the molecular solvent SPC water, and the polyion with the helical charge distribution. The results are the averages over 15 configurations, selected every 5 X lo5 trial steps, from run D. Upper figure: ( 0 )Total electrostatic potential. The solid curve, included for comparison, is for the axially charged polyion and the continuum solvent (c = 66). It is the same as the solid curve in Figure 9 (upper). Lower figure: breakdown of contributions to the total electrostatic potential. ( 0 )Total electrostatic potential.
counterions in run E, getting stuck in a local potential energy minimum just outside the glassy wall. It did not happen in run D where all 15 counterions clustered near the polyion. The singlet distribution function of SPC water for this model is shown in Figure 11. The narrow high first peak seen in Figure 3 for the uniformly charged polyion model is now absent. In fact this distribution function (Figure 11) somewhat resembles the one shown in Figure 4, for SPC water around an uncharged polyion. Clearly, for the helical necklace model of the polyion, there is extensive screening of the negative field of the fixed ions by the counterions, so that now the SPC particles feel a much weaker
1 l.0(rwp(12.0A
12.0(rwp(l3.OA
electrostatic field from the polyion than was previously the case. This is borne out as well by the angular distribution functions P,,(B) and PoH(B)shown in Figure 12. A comparison of Figures 12 and 6 (note especially the different scales for the vertical axes in the upper panel) reveals the much milder degree of dipolar solvent polarization with the necklace model, relative to the uniformly charged model. The large peak in POH(6)a t 8 < rwp < 9 and 170 < B < 180 is indicative of hydrogen-bonding of the SPC solvent with a negative charge on the polyion surface. As expected, this decays rapidly with distance. Finally, we display the electrostatic potential throughout the simulation cell for the necklace model in Figure 13. Our result is noisy because it is based on only 15 configurations. Despite the uncertainty, it is seen that the electrostatic potential is negative near the polyion surface, which leads to the counterions clustering near it. The solvent no longer contributes as large a positive electrostatic potential as it did with the axially charged polyion (Figure 9, bottom). In summary, we have found that the use of a uniform linear charge density along a central axis (which is electrostatically equivalent to a uniformly charged cylinder) as the representation of the charge distribution of a polyion leads to an unrealistically high degree of solvent polarization, when SPC water is used for the solvent. This leads to a counterion distribution wherein the counterions avoid regions close to the polyelectrolyte. When the polyion's charge distribution was made more realistic, by modeling it as a helical necklace of charges at the surface of a cylinder, the degree of solvent polarization was much reduced, and the counterions now formed a compact double layer near the surface of the polyion. Clearly, the form taken by the counterion distribution function, near a polyion, is a very delicately balanced function, which can change dramatically with changes in the model for the polyion. This should be borne in mind in future studies when polyion models are considered.
Acknowledgment. We thank Peter Rossky, Martin Neumann, and Tony Haymett, whose comments, taken collectively, had a significant and salutary effect on the final form of this work. We also thank Tom Record for introducing us to this subject and the Natural Sciences and Engineering Research Council of Canada for financial support. Appendix Potential Energy Expressions. The long-range nature of the Coulombic part of the potential energy makes it necessary to include these contributions from distances beyond the axial extremities of the central cell. It is in fact necessary to include these contributions, albeit approximately, from distances extending to infinity above and below the cell. For this reason we split up the
1930 The Journal of Physical Chemistry, Vol. 96, No. 4, 1992
potential energy expressions into “internal” and “external” contributions. Internal contributions are those due to interactions at distances within the central simulation cell; external contributions are from interactions over distances beyond the vertical extremities of the cell. (a) Polyelectrolyte with Uniform Linear Charge Distribution in a Dielectric Continuum Solvent. The total potential energy of an ion i for a particular configuration was obtained from
u, = ullnt+ ulext
Gordon and Goldman
tb The term uipEL,ext(rip) in eq A3 was obtained in a similar way. As in other work20 UiDEL,ext =
(AI)
where uIlnt =
(A9)
N
E ( ~ l l L J ( r l+J )~ ~ ~ ~+ ulpLJ(rIp) ~ ( r + ~ ~~l p )E L) T r l p )rIp t < Ro
J=l
Replacing the infinite limits -m and +- in eq A9 by -H and +H, respectively, such that H >> Ro and H >> LJ2, and doing the integrations gives
J#I
--
(A2) m,
ripIRo
and
+
Ulext = ulpEL~ext(rlp)uIiext(rlp), rlp< Ro
-- m,
(A3)
rlpIRo
The superscripts EL, LJ, int, and ext stand for electrostatic, Lennard-Jones, internal, and external, respectively, and the subscripts ii and ip mean ion-ion and ion- lyion, respectively. Ro is the radius of the simulation cell (25 r F i g u r e l), N is the total number of counterions in the cell (1S), rlpis the radial distance between ion i and the polyion axis, and rlJis the center-to-center distance of separation between counterions i and j . The smooth cylindrical outer wall of the simulation cell is created by the condition U,= m for rlpI&. The Lennard-Jones and electrostatic contributions to VIntdue to interactions with other counterions were obtained from
The first term on the right-hand side of eq A10 is divergent, but, as will be shown, cancels with a similar divergent term that arises in uiTt. Since the form of the counterion distribution function above and below the central cell is not known, an approximation is needed for the term uiiat(rip)in eq A3. As in other studies on the electrical double and on polyelectrolytes:0 where this problem arose, we approximated the external counterion distribution by setting it equal to the internal counterion distribution, which is generated and periodically updated during the simulation. Specifically, we used
where pext(r‘) is the average charge density due to counterions above and below the cell. Going over to cylindrical coordinates, replacing the infiite limits in the axial direction by f H , and using the electroneutrality condition in the form
eb = 2 r L b pext(r’)r’dr‘
and
uiPL = qiqj/erij
(A51
where qiLJand uikJare the Lennard-Jones parameters (Table I), qi = zle,where zi is the ionic valence (zi = zj = l), e is the electronic charge, and e is the static dielectric constant of the continuum solvent. From eq A2, a counterion interacts with the segment of the polyion inside the simulation cell, through a Lennard-Jones-like term, and an electrostatic term. The Lennard-Jones-like term was obtained by uipLJ
= 4eipLJ[
The values used for
ti:J,
(
I):(
(A6)
?$)12-
uitJ,in eq A4 and
tipLJand
gives 2zie2 cb
uiFt = -log, (2HJ- 2z1e e
LRo pext(r’)r’dr’
X
(A131 The divergent first term in eq A13 cancels the divergent first term in eq A10, when the two expressions are combined to create a working equation for the potential. Note also that, after eliminating the divergent term in eq A10, the total ion-polyion potential is found, by adding eqs A8 and A10, to be
usLJ in eq
A6 were those which made our system behave like a polyelectrolyte
system studied e1se~here.I~ Specifically, we wanted the minimum in the counterion-polyion potential (Coulombic plus non-Coulombic parts) to occur at rip = 10 A. This was achieved by using the parameters entered in Table I. The ion-polyion electrostatic potential energy of interaction inside the simulations cell is obtained by integrating the Coulombic potential between the point charge, qi, on the ion and the line charge of density -e/b over the axial length, L, of the simulation cell:
Doing the integration gives
The interactions included in this calculation are summarized by the entries at the top of Table 11. (b) Polyelectrolyte with Uniform Axial Charge Distribution in SPC Water. Here we change the model from the one described above by replacing the dielectric continuum inside the simulation cell by SPC water molecules and by replacing the smooth outer cell wall at & by the glassy wall described previously. The purpose of the glassy wall was to reduce wall-induced solvent structure that smooth walls next to molecular solvents are known to induce.37 This is discussed extensively in ref 37. We again approximate the external counterion distribution by an average charge density obtained from the charge density inside the cell, but there now arises the additional complication of what to do about the solvent above and below the cell. The traditional way of handling this
The Journal of Physical Chemistry, Vol. 96, No. 4, 1992 1931
Counterion and Solvent Distribution Functions is to do an Ewald sum,33 but this would have been too timeconsuming a calculation, because of the large number of solvent particles in our central cell (1454 mobile SPC's and 1465 frozen SPC's in the glassy wall). We handled this problem by doing separate calculations using two different approximations to model the effect of the external solvent. We call these approximations, which were designed to bracket the actual solvent effect, the "overscreenedn and Yunderscreenednapproximations. In the overscreened approximation, Coulombic interactions beyond the vertical extremities of the cell are screened by a factor of l/e (with t = 66). Doing this clearly overscreens these external interactions because they are simultaneously also being screened by the polarized SPC molecules inside the cell. Similarly, our underscreened approximation involves no screening of distant interactions, beyond that due to polarized SPC particles inside the central cell. We reasoned that because these approximations bracket the actual degree of screening of distant Coulombic interactions, the results of the two calculations should bracket the correct answers. The interactions included in these calculations are summarized by the entries in the middle of Table 11. The expressions we used for them are given below. Consider the internal interactions first. A counterion interacts with other counterions, the polyion, the SPC waters (both mobile and those in the glassy wall). The counterion-counterion, and counterion-polyion potential energy expressions are given by eqs A4, A5, A6, and A8. The dielectric constant e in eqs A5 and A8 is now set equal to unity and the Lennard-Jones parameters used in eqs A4 and A6 are given in Table I. The counterion-water interaction between an ion i and water molecule j (frozen or mobile) was taken to be
ujw =
#jWLJ
+ #jwEL
where
(
uiwLJ= 4tiwLJ[ $)Iz
-
I)$(
(A19
(A16)
The values of and u,~" were taken to be those of SPC water (Table I). Equation A17 gives the Coulombic part of the interaction energy of an ion with an SPC water molecule. q," is the ath partial charge at r? on water moleculej. The water-water interaction energy was obtained from the SPC pair potential function (ref 31 and Table I), with a spherical cutoff at 2 . 9 4 0 ~ ~ . Interactions between frozen SPC particles within the glassy wall or between frozen SPC particles and the polyion were not needed because these pair energies never changed during the simulation. The water-polyion interaction energy inside the simulation cell was obtained using int = uwpLJ + UwpEL,int ('418) UWP
were approximated by an average continuous charge distribution. The Coulombic interaction between this charge distribution and an ion (at z = 0) in the cell was screened by a dielectric constant of 1 or 66, depending on whether our simulation was for the underscreened or overscreened approximation. Thus we used eq A13 for this calculation, with e taken as 1 for the underscreened approximation and 66 for the overscreened approximation. The external contribution to the ion-polyion interaction energy was obtained by eq A10 with e taken as 1 and 66 for the underscreened and overscreened approximations, respectively. The external contribution to the water-polyion potential was obtained by
and the external contribution to the water-counterion potential was obtained by
[(tr+
? -
(rwpa)2+ r r Z- 2rwpur' cos d ] ' / 2d}d (A22)
where e was 1 or 66 depending on whether the underscreened or overscreened approximation was being used. Equations A21 and A22 are obtained in the same way as eqs A10 and A13. (c) Polyelectrolyte with Discrete Helical Charge Distribution in SPC Water. The solvent and counterions in these calculations were the same as in part b (above). Also, the polyion segments above and below the cell were as in part b. Each was a smooth soft cylinder extending to plus or minus infinity with a uniform axial charge distribution. The model of the polyion used here differed from the one in part b only in that here the charge distribution on the polyion segment inside the simulation cell was a helical necklace of discrete univalent negative charges, as opposed to a uniform axial charge density. Consequently, all external interactions were obtained exactly as in part b. Also, in these calculations, we used only the overscreened approximation, so that the value o f t in eqs A10, A13, A21, and A22 was always 66. The only potential energy expressions that were modified in these calculations were those for the internal contributions to the counterion-polyion and SPC-polyion interactions. To maintain the same overall charge density on the polyion as in the models described in parts a and b, we had to distribute 15 univalent negative ions on the surface of a cylinder of height 25.5 A. We also wanted the polyion radius to be 10 8, and the spacing to be sufficient to allow for possible hydration shells to exist around each of these ions. These requirements were met by placing 15 univalent negative ions on a helix, 2 turns of which are contained by our simulation cell. The equations for the coordinates of these fixed ions are
zi = 0.85 + 1.7(i - 1); i = 1, 2, ..., 15
where xi = 10 sin
(A20) Here, rwpis the radial distance between the polyion axis and is the radial distance between the center of the SPC water; rWp" the partial charge qwaon the water molecule and the polyion axis. Equation A20 was obtained in the same way as eq A8. The external interactions we included are given by the entries in the third column, in the middle of Table 11. As with the continuum solvent calculations, ions outside the simulation cell
$i,
yi = 10 cos (jJi
(A23)
In these equations (xi, yi, zi) are the space-fixed Cartesian coordinates, in A, of the 15 charges. This charge distribution has each negative ion centered on a cylindrical surface 10 A from the polyion or cell axis. The center-to-center distance between ions is 8.31 A. As already mentioned, the only potential energy expressions that change are those that describe interactions inside the cell for the polyion-water and polyion-counterion interactions. Each of the 15 negative charges on the polyion surface was contained at the center of a Lennard-Jones sphere whose u and e parameters were those of SPC water (Table I). These charged Lennard-Jones spheres crudely mimic charged functional groups on a polyion. Consequently, eq A2, was replaced here by
J. Phys. Chem. 1992, 96, 1932-1938
1932
Jones and electrostatic terms in eq A15 are of the form given in eqs A4 and AS, with the dielectric constant in eq A5 taken as J#i
('424)
unity. Similarly eqs A19 and A20, for the internal contribution to the water-polyion potential for these calculations, were replaced by
u,",LJ=
where I5
uipint(ri)=
C (uimLJ(lri- rml) + uimEL(lri- rml)]
m= I
('425)
The sum in eq A25 is over the 15 fixed negative ions on the polyion segment inside the cell. Ir, - rml is the center-to-center distance of separation between counterion i and fixed polyion charge m. Notice that we now have two types of Lennard-Jones terms in (A24) and (A25): uipLJ(rip) and uimu(lri- rml). The first, eq A6, prevents counterions from penetrating too far into the body of the polyion cylinder, and the second prevents them from getting too close to the negative charges on the polyion. The Lennard-
and
respectively.
Experimental and Theoretical Studles of the Equation of State of Liquid 2-Butyne from 248 to 293 K and Pressures up to 104 MPa V. Garcia Baonza,* M. Caceres AIonso, and J. N l e z Delgado Departamento de Qdmica F k c a , Facultad de Ciencias Qujmicas, Universidad Complutense, 28040 Madrid, Spain (Received: July 22, 1991)
The expansion method has been used to measure the density of liquid 2-butyne, from 248 to 293 K and pressures up to 104 MPa. The solid-liquid coexistence curve up to 100 MPa has also been determined. The experimental results have been fitted to the Strobridge equation of state, which has been used to estimate several properties of compressed liquid. The abilities of other semitheoretical and theoretical equations of state to predict the density and thermodynamic properties have been examined and discussed.
Introduction Accurate PVT data of simple liquids over a wide range of densities are known to be suitable for testing semitheoretical and theoretical equations of state and very useful for improving intermolecular potential functions.'" In the last years, equations of state based on statistical mechanics have been proposed in order to attempt a molecular understanding of the properties of a fluid. Some of the more interesting equations of state can be obtained from perturbation theory, and we will concentrate our attention on a semitheoretical one, the D e i t e ~ sequation ~,~ of state, which has been used to predict thermodynamic properties of 2-butyne over the whole density range, including liquid and vapor phases up to the critical point. Another theoretical equation of state, developed by Haar and Kohler,6 has been tested in terms of our experimental density results and has been used to estimate thermodynamic properties at low densities, such as second virial coefficients, from highpressure results. In this work we present the experimental study of the equation of state of liquid 2-butyne a t pressures from saturation to about 104 MPa and temperatures from 248 to 293 K. This substance is very interesting from the theoretical point of view because of its high degree of linearity and rigidity that allows one to compare (1) Calado, J. C. G.;Clancy, P.; Heintz, A,; Streett, W. B. J . Chem. Eng. Data 1982, 27, 376. (2) Rubio, R. G.; Calado, J. C. G.; Streett, W. B. In Equations of State: Theory and Applications; Chao, K.C., Robinson, R. L., Eds.; ACS Symposium Series 300; American Chemical Society: Washington, DC, 1986. (3) Powles, J. G.; Evans, W. A.; McGrath, E.; Gubbins. K. E.: Murrad. S . Mol. Phys. 1979, 38, 893. (4) Deiters, U. Chem. Eng. Sci. 1981, 36, 1139. (5) Deiters, U. Chem. Eng. Sci. 1981, 36, 1147. (6) Kohler, F.; Haar, L. J . Chem. Phys. 1981, 75, 3 8 8 .
it with theoretical calculations. In addition, this work brings new and accurate experimental data at high liquid densities, up to the melting line, which are very scarce in the literature. The results consist of 111 PpT points measured along seven isotherms, which covers a relatively wide region not previously studied, it is important to remark that the 2-butyne remains liquid only from 241 to 300 K at atmospheric pressure. The experimental densities have been correlated in terms of the Strobridge' equation of state, which has been used to calculate the isothermal compressibility and thermal expansion coefficient at regular intervals of pressure and temperature. The experimental study has been completed with the measurement of the solid-liquid coexistence curve for pressures below 100 MPa which has been used, together with vapor pressure and critical data, to construct the phase diagram and estimate the triple point of this substance. The second virial coeficient as a function of temperature has been calculated from calorimetric and vapor pressure data, and the results have been compared with those predicted by semitheoretical and empirical expressions as well as those derived from the theoretical equations of state indicated above. Special attention has been given to the behavior of this property at low temperatures.
Experimental Section The results reported here have been obtained with an expansion technique. A detailed description of the method and apparatus used in this work has been published.'.8 Temperatures have been measured with an accuracy of 0.01 K,using a Leeds and Northrup calibrated platinum resistance thermometer, and were referred to the international temperature (7) Strobridge, T. R. NBS Tech. Note (US.)1962, No. 129. (8) Baonza, V. G.; Niifiez, J.; Cieeres, M. J . Chem. Thermodyn. 1989,21, 231.
0022-3654192 12096- 1932%03.00/0 , 0 1992 American Chemical Society I
,
