5714
J . Phys. Chem. 1987, 91, 5714-5718
Brownian Dynamics Simulation of Counterion Dynamics In Cylindrical Polyelectrolyte Solutions Lars Guldbrand* and Lars Nordenskiold Division of Physical Chemistry, Arrhenius Laboratory, University of Stockholm, S-106 91 Stockholm, Sweden (Received: January 15, 1987)
We report an investigation of cylindrical polyelectrolyte systems using the Brownian dynamics computer simulation method. The counterion dynamics is studied, and the results are compared with the predictions of m a n field theory. We find a significant but highly model dependent retardation of the counterion motion.
Introduction The theoretical description of charged macromolecules in solution is often made using the primitive model, the cell model, and the mean field approximation.' The mean potential in the vicinity of the macroion is then given by the solution of the Poisson-Boltzmann (PB) equation in the appropriate symmetry. In studies of ion dynamics, the basis is often the combination of this PB potential with the Smoluchowski equation.2 Much work has recently been performed concerning the applicability of the PB equation in different polyelectrolyte systems.- In some cases the neglect of the small-ion correlation has been shown to result in significant errors.lOql' A first step to improve on the mean field Smoluchowski description can then be to use a potential of mean force calculated by, for instance, computer simulation methods and which takes small-ion correlations into account.I2 A more simple procedure would be to calculate an approximate potential of mean force from a solution of the modified Poisson-Boltzmann equation. For making more direct studies of dynamical properties, the molecular dynamics (MD) computer simulation method14 is an important tool. However, a good description of the system, including the solvent, is needed, and the time steps used in the simulation must be small. This makes the simulation of slow dynamic processes very expensive in computer time. In many cases, some species in the system are of primary interest, while detailed information is not sought about the remainder. For instance, in the study of small-ion dynamics in polyelectrolyte systems, one might be prepared to stay within the primitive model and thus forsake detailed information about the solvent in exchange for being able to follow the motions of the counterions for long periods. In this case, the Brownian dynamics (BD)lS-I8 (1) Verwey, E. J. W.; Overbeek, J. Th. G. Theory of the Stobility of Lyophobic Colloids; Elsevier: Amsterdam, 1948. (2) Chandrasekar, S.Rev. Mod. Phys. 1945, 15, 1. (3) van Megen, M.; Snook, I. J. Chem. Phys. 1980, 73,4656. (4) Torrie, G. M.; Valleau, J. P. J . Chem. Phys. 1980, 73, 5807. (5) Jhsson, B.; Wennerstrom, H.; Halle, B. J . Phys. Chem. 1980, 84, 2179. (6) Bratko, D.; Vlachy, V. Chem. Phys. Lett. 1982, 90, 434. (7) LeBret, M.; Zimm, B. H.Biopolymers 1984, 23, 271. (8) Murthy, C. S.;Bacquet, R. J.; Rossky, P. J. J . Phys. Chem. 1985,89, 701. (9) Mills, P.;Andersson, C. F.; Record, M. T., Jr. J . Phys. Chem. 1985, 89, 3984. (10) Guldbrand, L.; Jonsson, B.; Wennerstrom, H.;Linse, P J . Chem. Phys. 1984, 80, 2221. (1 1) Kjellander, R.; Marcelja, S.Chem. Phys. Lett. 1984, J 12, 49. (12) Akesson, T.; Jonsson, B.; Halle, B.; Chan, D. Y. C. Mol. Phys. 1986, 57, 1105. (13) Bratko, D. Bioelectrochem. Bioenerg. 1984, 13, 459. (14) Kushick, J.; Berne, B. J. In Modern Theoretical Chemistry; Berne, B. J., Ed.; Plenum: New York, 1977; Vol. 6, Chapter 2. (15) Ermak, D. L. J . Chem. Phys. 1975, 62, 4189. (16) Turq,P.; Lantelme, F.; Friedman, H. L. J. Chem. Phys. 1977, 66, 3039. (1 7) van Gunsteren, W. F. H.; Berendsen, J. C.; Rullman, J. A. Mol. Phys. 1981, 44, 69.
0022-3654/87/2091-5714.$01.50/0
simulation method can be an interesting alternative to the conventional M D technique. In this work, we have studied the counterion dynamics in primitive models of cylindrical polyelectrolytes using the BD simulation technique. For comparison, we have also applied the Smoluchowski equation with both PB potentials and potentials of mean force calculated by simulation methods.
The Simulation Method The motion of a Brownian particle surrounded by solvent molecules can often be described by the strict Langevin equation2 (or equivalently by the Smoluchowski equationI2). For the limit of infinite dilution, and with no external field, we have
Here, p is the linear momentum, { is a frictional coefficient, and
17 is a random force due to collisions with the solvent molecules. The random force and the frictional coefficient are related through a form of the fluctuation-dissipation theorem:
The displacement of the particle during a time step At, which is long compared to l/{, can be written At AF(t,At) = A(t)mf
(3)
where a ( t ) is an effective random force arising from the many quasi-instantaneous interactions of the particle with the solvent molecules during At. In the limit of infinite dilution, this effective random force follows a Gaussian distribution of zero mean and with a second moment which is related to the diffusion constant Do by 6mkT!: 6(kq2 (lAl2) = -= (4) At DoAt
In a system of interacting Brownian particles, an external force field E can be introduced in the strict Langevin equation, and we can assume that the motion of a particle can be described by
The fluctuation-dissipation theorem in the presence of this external force can now be formulated'*
&-((R(O) R ( t ) ) )dt
+ &-e-rr([R(O) + P ( O ) ] E ( t ) )dt
(6)
The second term in eq 6 is often small. When this is the case, (18) Akesson, T.; Jonsson, B. Mol. Phys. 1985, 54, 369.
0 1987 American Chemical Society
The Journal of Physical Chemistry, Vol. 91, No. 22, 1987 5715
Cylindrical Polyelectrolyte Solutions the distribution of the effective random force can be approximated by the Gaussian distribution of the infinite dilution case. The interaction between the particles is then assumed to have no effect on the dynamic properties of the solvent. For a time step Af, which is much longer than l/{ but so short that F(f+At) = F ( f ) ,the displacement of a particle can be written in a simple w a ~ : ' ~ J * At Ar(t,At) = [ d ( f )+ p(f)](7) mT This relation can be used as the basic algorithm for a BD computer simulation of interacting Brownian particles. During each time step, the velocity correlation is completely lost, and the simulation can be carried out in configuration space.
Models In a BD simulation of a polyelectrolyte system described within the primitive model, the solvent averaged potential between the ionic species is obtained through the introduction of the relative dielectric permittivity into the expressions for the ionic interactions. The distribution of random forces can be generated by using measured self-diffusion constants in aqueous solution. We have used the primitive model and the cylindrical cell model symmetry to describe systems of mono- and divalent counterions in the presence of cylindrical macroions. Our main interest has been focused on the B-DNA double helix, but we have also studied monovalent counterions in the presence of poly(acry1ic acid) (PAA) at different degrees of neutralization. The macroions were pictured as hard cylinders of radius a = 10 8,for the DNA systems and a = 5 8,for the PAA systems, respectively. The counterions were modeled as point charges, and no co-ions were considered. To study the dependence of ionic size in an approximate way, a distance of closest approach, d, was used, so that the effective macroion radius was given by a + d. Three different models were used for the DNA charge distribution. In the first, model A, the macroion is considered to have a uniform surface charge of u = -0.15 C/mZ. In model B, the macromolecular charge is distributed as discrete elementary charges equidistantly in line on the DNA axis, while model C is an attempt at a more realistic charge distribution within the macroion cylinder. The unit charges are placed at the approximate locations of the phosphorus atoms in two helices of radius 9 8,. Each helix has a step per charge of 3.4 8, in the axial direction and 36' in the azimuthal plane, and the two helices are offset by 0.825 8, and 206°.19~20 For the PAA systems, only a discrete macroion charge distribution was used. In this description, model D, the unit charges are placed in a single helix of radius 4.0 8,with a step of 2.5 8, in the axial direction and 90' in the azimuthal plane. The simulations in the DNA systems were performed at a temperature of 298.15 K, and the solvent water was characterized by a relative dielectric permittivity of 78.3. The monovalent counterion considered was Li+ and the divalent Ca2+. This choice was determined by the fact that the self-diffusion of these ions in DNA solutions has been studied experimentally, albeit under different solution conditions than the present simulations.2' The infinite dilution diffusion constants were taken to be 1.07 X and 0.79 X m2/s, respectively.21q22 For the PAA systems, we used a dielectric permittivity of 77.94, which corresponds to D 2 0 at 298.15 K. The counterion considered is again Li', but the self-diffusion constant at infinite dilution in DzO is 0.86 X 1 0 - ~m2/s.23 For the DNA systems, the simulations were made in cylindrical simulation cells with a radius, r, of 50.0 8, and a length, I, of 136 ~~
~
(19) Clementi, E. In Lecture Notes in Chemistry; Springer-Verlag: Berlin, 1980; Vol. 19 and references cited therein. (20) Guldbrand, L.; Nilsson, L. G.; Nordenskibld, L. J. Chem. Phys. 1986, 85, 6686. (21) Nilsson, L. G.; Nordenskiold, L.; Stilbs, P.; Braunlin, W. H. J . Phys. Chem. 1985.89, 3385. (22) Magdalenat, H.; Turq, P.; Chemla, M. Biopolymers 1974, 13, 1535. (23) Rymden, R.; Stilbs, P. J . Phys. Chem. 1985, 89,2425.
TABLE I: Cell Averaged Axial and Azimuthal Counterion Diffusion Quotients from the Simulations of the Different Charge Distribution Models of DNA" q 1 1 1 1 1
model
d, A
A
0.0 0.0 0.0 3.0 3.0
0.81 0.81 0.53 0.91 0.89
f 0.02 f 0.06 f 0.11 f 0.02 f 0.03
D,lDO 0.84 f 0.02 0.80 z t 0.04 0.53 f 0.07 0.92 f 0.02 0.92 f 0.02
2 2 2 2
A
0.0 0.0 3.0 3.0
0.28 0.46 0.41 0.46
f 0.20 f 0.15 f 0.15 f 0.20
0.41 0.59 0.61 0.57
" q is the
B C A C
B A C
DJDO
f 0.19 f 0.13 f 0.08 f 0.09
counterion valency, and d is the distance of closest ap-
proach. TABLE II: Cell Averaged Axial and Azimuthal Counterion Diffusion Quotients from Simulations of the Model D of Poly(acrylic acid) with Different Degrees of Neutralization deg of q model neutr DJDO DdDO ~~
1 1 1
D D D
100% 75% 50%
0.88 f 0.05 0.84 f 0.02 0.90 f 0.02
~
0.93 f 0.02 0.90 f 0.03 0.94 i 0.03
8,. For the PAA systems, the values r = 46 8, and I = 200 8, were used. These sizes correspond to 80 unit charges on the macroion. For PAA, the cell volume corresponds to a monomer concentration of 0.100 M. The minimum image convention and periodic boundary conditions were used in the axial direction, while the cell surface was treated as a hard wall. As the electrostatic interaction is of very long range, it is desirable to include, somehow, the interactions with the charges outside the simulation cell in the axial direction. The charge distribution here is approximated by the corresponding mean field d i s t r i b u t i ~ ntaken ~ ~ ~from ~ the analytical solution of the Poisson-Boltzmann equation in cylindrical symmetry.z4 The time step used in the simulations was 0.2 ps. The simulations were generally started from random configurations weighted with the mean field distribution and were equilibrated for about 5-10 ns before the production runs. Production runs were 12.8 ns for the monovalent counterions in the DNA systems and 25.6 ns for the divalent, while the production runs of the PAA systems were 11.2 ns long. The applicability of the algorithm was tested by means of eq 6, and it was found that the second term was negligible for all the monovalent simulations. However, for D N A systems with divalent counterions, discrete macroion charges, and no distance of closest approach, the second correlation function integral term of eq 6 could not be ignored, and the simulation method broke down. This problem was not encountered when a distance of closest approach of 3.0 8, was used. Results and Discussion The equilibrium properties of the model systems are given by the mean distribution of the counterions. The simulated radial concentration (and potential of mean force) profiles are all rather similar to those of the corresponding PB calculations. The static correlation effects included in the simulations tend to raise the concentration at the macroion surface and lower it at the cell b ~ r d e r . ~This - ~ effect is, as can be expected, more pronounced for systems where the correlation is important, so that the systems with discrete wall charges are more affected than those with a uniform charge density and those with divalent more than those with monovalent counterions. Some examples of radial concentration profiles are displayed in Figure 1 . It can be seen that the choice of macroion charge model means relatively little for monovalent counterions and distance of closest approach d = 3 A. For divalent counterions, however, the discrete helical charge (24) FUOSS, R. M.; Katchalsky, A.; Lifson, S.Proc. Natl. Arad. Sei. U.S.A. 1951, 37, 579.
5716 The Journal of Physical Chemistry, Vol. 91, No. 22, 1987 4.0
x
40
0.3
a
c
3.0
K v
0.2
=
v
h
30
L h
L
L
v
2.0
0.1
.o
0.0
1
,I
a
h
v
"
I
I
l
VI
h
v
Guldbrand and Nordenskiold
t 0 L v
t-
20
10
0.0
0 20
30
r
5.0 4.0
50
40
(A)
0.10
')I h
I
I
h
t
M
0.08
c I
v
v
0.06
h
n L
3.0
0.04 y
L
t 0 L v
2.0
1
0.02
0.00
.o
0.0
t-
'. I
to
-._.
I
I
I
30
50
(A) Figure 1. Radial concentration profiles for (a) monovalent and (b) divalent counterions surrounding a DNA cylinder with a radius of 10 A and a distance of closest approach of 3 A. The solid lines represent the predictions of mean field theory, while the dashed and the dotted lines represent values from the model A and model C simulations, respectively. (The insets show the concentration for radii greater than 20 A in a larger scale.) r
distribution results in a concentration profile with more ions close to the DNA relative to that of the model A. Diffusion coefficients are usually evaluated by using the mean square displacement. However, hesson and JBnsson have derived an alternative expression involving the force correlation function (eq 40 of ref 18). This method was shown to give better numerical precision, and we have used it here to calculate the diffusion coefficients for the self-diffusion of the counterions. Tables I and I1 show the cell averaged diffusion constants D,and D, for the axial and the azimuthal motion of the counterions, respectively, for both the DNA and the PAA systems. As can be seen from the results of the monovalent DNA systems, the difference between the diffusion in models A and B is not great, while model C seems to give a more substantial retardation. The distance of closest approach is however an important factor, and when this is set to 3.0 A, the diffusion coefficients of the monovalent counterions are close to the infinite dilution value. For divalent counterions, both model A and model C give a considerable lowering of the diffusion coefficients. In Table 11, there seems to be a minimum in the diffusion coefficients for the 75% degree of neutralization. However, considering the estimated errors, this need not be a statistically significant result. For the same reasons, it is also difficult to be certain about the tendency of the azimuthal diffusion coefficient to be slightly larger than that of the axial motion in the same system. The dynamics in the radial direction were characterized by mean first passage times (MFPT)25with a cylindrical, completely reflecting, inner boundary at r, and a completely absorbing outer boundary at radius rb.26 The MFPT T(ro)is then the mean time (25) Akesson, T.; Jonsson, B. J . Phys. Chew. 1985, 89, 2401
20
30
40 r
50
(A)
Figure 2. Mean first passage times as a function of the location of the absorbing boundary for (a) monovalent and (b) divalent counterions in the DNA system with a 3-A distance of closest approach. The initial position r,, is placed at the effective radius of closest approach. The dashed line shows results from the simulations with the charge distribution of the model A, while the dotted line represents the results of model C. For comparison, the solid lines show values calculated by use of the PB-Smoluchowski model.
it takes for a particle initially located at the radius ro (ra 5 r, 5 rb) to become absorbed at radius rb. In the simulations, various MFPT's were calculated from r,, = ra = a d; Le., the inner boundary and the initial positions were both located at the surface of closest approach of the ions to the macroion. The outer boundary was placed at different radii within the cell. The production runs were rather short compared to the time scale expected for the MFPT's, and only the values for rb - ro less than about 5 A show no significant changes with increasing simulation length. For the larger values of rb, the MFPT's were instead calculated by using the Smoluchowski equation and the simulated, axially averaged, radial potentials of mean f 0 r ~ e . IFor ~ ~the ~~ smaller distances between r, and rb, this procedure resulted in values not significantly different from those directly calculated. This indicates that the effect of dynamical correlations is small in these systems. The MFPT's seem very sensitive to the form of the potential, and differences that show up as very small variations in the concentration profiles may result in large variations in T . The results for the DNA systems are shown in Figure 2 together with the predictions of the PB-Smoluchowski The differences between the two curves for the simulations of monovalent counterions using model A and C, respectively, correspond to the very small relative differences in the concentration profiles of Figure la, and they are probably not significant. For divalent counterions, on the other hand, it can be seen that the results are highly dependent on the model used for the ma-
+
(26) Chan, D. Y . C.; Halle, B. Biophys. J . 1984, 46, 387
Cylindrical Polyelectrolyte Solutions 1.0
*
The Journal of Physical Chemistry, Vol. 91, No. 22, 1987 5717 TABLE III: The Zero-Time Value of the Electric Field Gradient Correlation Function (V(O):V(O)) from Simulations of the Different Models of tbe DNA Charge Distribution for both Monovalent and Divalent Counterions, together with the Results of the Corresponding Mean Field Calculations" (V:V)l/Z,1016
x , kHz
V/m2
1
0.2 0'4
0.0
+
1
+ + +
0.2 0.4 0.6 0.8 1 . 0 DEGREE OF NEUTRALIZATION Figure 3. The Li+ macroscopic counterion self-diffusion quotient for an isotropic solution of poly(acry1ic acid) (PAA) as a function of the degree of neutralization. The dotted line shows the simulation results, while the solid line shows the mean field theory values. The experimental results (+) are taken from ref 23.
croion charge. The differences between the model A and C curves in Figure 2b are directly comparable to the results of Figure 1b for the concentration profiles. All the simulations show substantial differences relative to the mean field calculations. The macroscopic counterion self-diffusion in a polyelectrolyte system can be measured by using radioactive tracer methods2' or the FT N M R pulsed gradient self-diffusion method.28 Previously, theoretical calculations have been made for macroions like DNA and PAA, relating the macroscopic diffusion to the diffusive motion of the ions in the cylindrical cell model. For an isotropic solution of cylindrical and infinitly long macroions, we have the relation2'
-D= - -1 D!l Do 3 Do
+ -23 -DDoL
Here, D,, is the diffusion coefficient for motions parallel to the macroion axis, while D , describes the diffusion in a direction perpendicular to the macroion. In the mean field description the parallel diffusion is entirely unhindered, and DIlis equal to Do, while the second term is dependent on the mean potential profile. Equation 8 is then equivalent to the expression for the macroscopic diffusion derived by Yoshida using mean field theory.29 For the model D PAA systems, we have used our directly simulated axial diffusion coefficients for Dll. In order to calculate D,, we have used the Smoluchowski cylindrical cell model result given by eq 21-23 in ref 21. However, instead of the PB potential we used our simulated, axially averaged, potential of mean force in these numerical calculations. By using these values for Dlland D , in eq 8, we obtain an estimate of the diffusion quotient in the PAA system which can take the discrete nature of the macroion charge and parts of the correlation effects into account. This procedure is clearly approximate and not entirely consistent with the microscopic Brownian dynamics model, since eq 8 was derived from the phenomenological Fick's first law,21but the result should still be an improvement compared to the PB-Smoluchowski (PBS) model. The results are directly comparable to experimental measurements made by Rymden and S t i l b using ~ ~ ~ the FT N M R pulsed gradient self-diffusion method. The simulation results are shown in Figure 3, together with the corresponding experimental and mean field values. The simulations are found to give results significantly lower than the mean field calculations, but no qualitative agreement with the experiments is obtained. The lowering of the diffusion quotients as compared to the PBS results is due to both reduced axial (Dli)and transverse diffusion (0,). (27)Huizenga, J. R.;Grieger, P. F.; Wall, F. T. J . Am. Chem. SOC.1950, 72, 4228.
(28)Stilbs, P.;Lindman, B. J . Magn. Reson. 1982, 48, 132. (29) Yoshida, N.J . Chem. Phys. 1978, 69, 4867.
q
model
1 1 1 1 1
A
2 2 2 2
A B A
B C A
C
C
d, A 0.0
0.0 0.0 3.0 3.0
SIM 83 f 27 109 f 105 812f 76 78 f 67 72 f 45
PB 105 105 105 25 25
0.0 0.0 3.0 3.0
58 f 18 34 f 4 36 f 18 36 f 22
85 85 51 51
SIM -0.49 f 0.16 -0.65 f 0.60 -4.79 f 0.45 -0.47 f 0.40 -0.42 f 0.27
PB -0.614 -0.614 -0.614 -0.143 -0.143
30.4 f 17.0 f 17.9 f 17.7 f
44.065 44.065 8.128 8.128
8.9 2.7 8.9 11.8
Values of the quadrupolar coupling constants computed from the field gradients are also given: x = (1 y.)eQ/h[(2/3)(V:V)]'/2. The commonly used polarization factor in this equation has been shown to be questionable and has therefore been omitted.@ The quadrupole and 0.2 X moments Q of 'Lit and 43Ca2+are taken to be -0.04 X m2, respectively, while the Sternheimer antishielding factors (1 y.) are taken to be 0.75 for Li+ and 13.1 for Ca2t.39
+
+
It should be mentioned that due to the coiled nature of vinylic polyions (particularly at low charge densities) like PAA, eq 8 overestimates the contribution from Dllto D. This is perhaps the most likely reason for the discrepancy between the theoretical and . experimental values. The measurement of the nuclear magnetic relaxation time offers an important way of investigating the surroundings of ions possessing a nuclear quadrupole moment.30 This relaxation rate arises from the interaction of the quadrupole with the electric field gradient at the ion nucleus. For ions in polyelectrolyte solutions, one could distinguish between two contributions to the electric field gradient: one arising from the polar solvent molecules in the neighborhood of the ion and the other due to the electric charges of all the ionic species in the ~ y s t e m . ~ ~The ' ~ first contribution has been studied by, for instance, Monte Carlo and molecular dynamic simulations and has been found to be an important relaxation mechanism for quadrupolar ions in aqueous solutions.33 The second contribution can be estimated in a continuum model, and this has been done by using mean field theory.34 This yields rather small field gradients compared to those due to the solvent. Nevertheless, attempts have been made to interpret measurements of relaxation rates using only the field gradient our main calculated from the PB d i s t r i b ~ t i o n . ~ ~ -Although ~' interest has been the counterion diffusion, we have also evaluated the zero-time value of the electric field gradient in the laboratory frame (VL(0):VL*(O)).The results are given in Table 111, together with values of the quadrupolar coupling constants, x, calculated for 7Li+and 43Ca2+,respectively. The coupling constant from the solvent contribution to the field gradient at a Li+ nucleus in ~ simulation ~ aqueous solution is of the order of -30 ~ H ZIn.the of the model C system with no extra distance of closest approach, (30)Abragam, A. The Principles of Nuclear Magnetism; Clarendon: Oxford, 1961. (31)Hynes, J. T.;Wolynes, P. G. J . Chem. Phys. 1981, 75, 395. (32)Halle, B.; Wennerstrom, H.; Piculell, L. J . Phys. €hem. 1984, 88, 2482. ( 3 3 ) EngstrBm, S.; Jonsson, B.; Impey, R. W. J . Chem. Phys. 1984, 80,
5481. (34)Wennerstrom, H.; Lindman, B.; Lindblom, G.; Tiddy, G. J. T. J. Chem. SOC.,Faraday Trans. 1 1979, 75, 663. (35)van der Klink, J. J.; Zuiderweg, L. H.; Leyte, J. C. J . Chem. Phys. 1974, 60, 2391. (36)Delville, A.; Gilboa, H.; Laszlo, P. J . Chem. Phys. 1982, 77, 2045. (37)Delville, A,; Laszlo, P. Biophys. Chem. 1983, 17, 119. (38)Engstrom, S.;Jonsson, B. Mol. Phys. 1981, 43, 1235. (39)Lucken, E. A. C. Nuclear Quadrupole Coupling Constants; Academic: London, 1969. (40)Engstrom, S.;Jonsson, B.; Jbnsson, B.J . Magn. Reson. 1982,50, 1.
5718
J . Phys. Chem. 1987, 91, 5718-5725
we find a field gradient far larger than the mean field estimate but still considerably smaller than that arising from the ionsolvent interaction. The spin relaxation is roughly determined by a product of the interaction and the time scale of its fluctuation. Since the solvent contribution has a characteristic time scale of picoseconds38 while the charge contribution fluctuates on a time scale of nanoseconds, it is not clear which will dominate the relaxation behavior. Conclusions The simulations of the DNA system show that, for monovalent counterions, there is a small but statistically significant retardation of both the radial as well as the axial diffusion, due to correlation effects. With a discrete double-helical charge model, the effect is larger but quite sensitive to the distance of closest approach, which is a somewhat arbitrary parameter. Experimentally measured Li+ self-diffusion coefficients in DNA solutions have been found to be slightly smaller than the predictions of the PB-Smoluchowski model.21 The present study then indicates that continuum model electrostatics could be able to give a reasonable description of the dynamics for monovalent counterions in this system and that the mean field treatment introduces only a rel-
atively small error.2o For divalent ions, the retardation is more significant. The results for the axial diffusion indicate that the low self-diffusion coefficients observed in solutions of NaDNA titrated with small amounts of Ca2+ might be described within the primitive model. For PAA, the calculated effective macroscopic diffusion coefficients are not low enough to be comparable to the experimental results. This may partly be due to the approximative method of evaluation. Perhaps more important, the model itself, picturing PAA as a rigid, impenetrable cylinder with a fixed helical charge distribution, may be erroneous in the present context, since PAA is a more linear and rather flexible molecule, particularly at low degrees of neutralization. In other respects, the continuum model is too simple. This is evidently mainly in the calculations of the electric field gradient, where the contribution arising from a molecular description of the solvent has to be taken into account.
Acknowledgment. We are grateful to Bo Jonsson, Bertil Halle, and Torbjorn Akesson for valuable discussions. This work has been supported by the Swedish Natural Science Research Council. Registry No. PAA, 9003-01-4; Li, 7439-93-2,
A Theory of Association in Fluids Composed of Linear Bifunctional Molecules Capable of Hydrogen Bonding Richard E. Boehm* and Daniel E. Martire Department of Chemistry, Georgetown University, Washington, DC 20057 (Received: April 10, 1987)
A mean-field lattice model theory of associated fluids comprised of linear bifunctional molecules which are assumed capable of participation in a maximum of two hydrogen bonds is formdated. Structurally the fluid is assumed to consist of an isotropic, heterogeneous mixture of either rodlike and/or chainlike associated clusters of various lengths and in general tortuosities which are continually disrupted and regenerated via rotational diffusion of the constituent molecules. A generalization of the DiMarzio lattice statistics developed by Herzfeld is utilized to formulate the configurational free energy of the fluid. Nearest-neighbor associative interactions which correspond to strong and weaker hydrogen bonds as well as nearest-neighbor dispersion interactions represent the internal energy contributions to the free energy. An equation of state is derived which is a generalization of the Ising lattice fluid equation of state obtained by Sanchez and Lacombe in that the average cluster size is determined by the average degree of association which is temperature and density dependent and determined from free energy minimization. The equation of state generates critical point compressibility factors which are less than 0.3 which agrees with experiment. When applied to liquid HF, the equation of state generates predictions of the critical compressibility factor, critical pressure, and density in excellent agreement with experiment. Explicit expressions for the heat capacity, isothermal compressibility, and thermal expansion coefficient are derived and results are calculated for a condensed and supercritical fluid chosen to simulate HF. The calculated heat capacity and its temperature dependence correlate well with the experimental results obtained for liquid HF. Both calculated and experimental heat capacities of liquid HF exhibit a monotonic increase with increasing temperature which is consistent with the description that heat absorption generates transitions from strong to weaker hydrogen bonds.
I. Introduction The present article describes the formulation of a computationally tractable statistical thermodynamic theory of simple associated fluids which are composed of bifunctional molecules which are capable of participation in the formation of a maximum of two hydrogen bonds with nearest-neighbor molecules. Condensed, supercritical, and vapor fluid phases of HF, HCN, and HC1 represent potential prototype systems. Structurally the fluid is assumed to consist of an isotropic, heterogeneous mixture of hydrogen-bonded chains of various lengths and in general tortuosities which are continually disrupted and regenerated via rotational diffusion of the individual constituent molecules. It is also assumed that the probability of hydrogen-bond formation between any nearest-neighbor molecular pair along a chain is completely independent of its position along the chain. The probability of hydrogen-bond formation is determined by free 0022-365418712091-5718$01.50/0
energy minimization and is a function of the temperature and density of the fluid and the strength of the association. The computer simulation study performed by Klein and McDonald' on condensed HF supports the picture that the liquid is characterized by a significant degree of hydrogen-bonded chain polymerization without ring formation. A mean-field approach based on a lattice model is utilized to formulate the free energy of the associated fluid. In particular we determine the number of accessible configurations and hence configurational entropy by adopting a generalization of the lattice statistics introduced by DiMarzio2s3and extended by Herzfeld4 (1) Klein, M. L.; McDonald, J. R. J . Chem. Phys. 1979, 71, 298. (2) DiMarzio, E. A. J . Chem. Phys. 1961, 35, 658. (3) Peterson, H. T.; Martire, D. E.; Cotter, M. A. J . Chem. Phys. 1974, 61, 3547. (4) Herzfeld, J. J . Chem. Phys. 1982, 76, 4185.
0 1987 American Chemical Society
