Charge

9600

Langmuir 2003, 19, 9600-9612

Overcharging in Macroions. Effects of Macroion Geometry/ Charge Distribution Arup K. Mukherjee,† K. S. Schmitz,‡ and L. B. Bhuiyan*,† Laboratory of Theoretical Physics, Department of Physics, University of Puerto Rico, Box 23343 UPR Station, San Juan, Puerto Rico 00931-3343, and Department of Chemistry, University of MissourisKansas City, Kansas City, Missouri 64110 Received April 15, 2003. In Final Form: August 15, 2003 A random sampling technique is employed to obtain the ground-state configuration energy of an overcharged macroion in that the magnitude of the total counterion charge exceeds the magnitude of the macroion charge. The technique is based on minimizing the total electrostatic energy of the macroion and counterions in salt-free conditions. The algorithm easily reproduces the ground-state energies of an overcharged isolated spherical macroion obtained from Molecular Dynamics simulation in earlier published literature (Messina, R., Holm, C.; Kremer, K. Phys. Rev. E 2001, 64, 021405) for low to intermediate values of the macroion charge, while at higher macroion charge the calculated energies lie below those of the earlier results. The calculations are extended to nonspherical macroion geometries such as spherocylinders and spheroids, and substantial effects on the energetics due to the geometries are noted. A detailed study of counterion distributions is also done for macroions of different charges and geometries.

1. Introduction Charge inversion of a macroion in an aqueous solution can occur when oppositely charged and usually simple small particles accumulate in the vicinity of the macroion surface due to electrostatic interaction in such an amount that the magnitude of the total counterion counterion charge becomes larger than the magnitude of the charge of the macroion itself. An early experimental confirmation of charge reversal was seen in the work of Strauss, Gershfeld, and Spiera.1 In these studies the electrophoretic mobility of “polysoap” derived from poly-4-vinylpyridine was reported to change direction in going from low concentrations of KBr (positively charged) to high concentration of KBr (negatively charged). This charge reversal was attributed to the excess binding of Br-. The phenomenon of charge reversal, or “overcharging”, of macroions has recently been suggested2-13 as a possible mechanism to explain experimental evidence that “likecharged” macroions appear to attract each other under certain conditions.14-17 † ‡

University of Puerto Rico. University of MissourisKansas City.

(1) Strauss, U. P.; Gershfeld, N. L.; Spiera, H. J. Am. Chem. Soc. 1954, 76, 5911. (2) Perel, V.; Shklovskii, B. Physica A 1999, 274, 446. (3) Shklovskii, B. I. Phys. Rev. E 1999, 60, 5802. (4) Nguyen, T. T.; Grosberg, A. Y.; Shklovskii, B. I. Phys. Rev. Lett. 2000, 85, 1568. (5) Nguyen, T. T.; Grosberg, A. Y.; Shklovskii, B. I. J. Chem. Phys. 2000, 113, 1110. (6) Mateescu, E. M.; Jeppesen, C.; Pincus, P. Europhys. Lett. 1999, 46, 493. (7) Joanny, J. F. Eur. Phys. J. B 1999, 9, 117. (8) Gurovitch, E.; Sens, P. Phys. Rev. Lett. 1999, 82, 339. (9) Deserno, C.; Holm, C.; May, S. Macromolecules 2000, 33, 199. (10) Messina, R.; Holm, C.; Kremer, K. Phys. Rev. Lett. 2000, 85, 872. (11) Messina, R.; Holm, C.; Kremer, K. Europhys. Lett. 2000, 51, 461. (12) Messina, R.; Holm, C.; Kremer, K. Eur. Phys. J. E 2001, 4, 363. (13) Messina, R.; Holm, C.; Kremer, K. Phys. Rev. E 2001, 64, 021405. (14) Ise, N.; Okubo, T.; Ito, K.; Dosho, S.; Sogami, I. Langmuir 1985, 1, 176. (15) Ito, K.; Ise, N.; Okubo, T. J. Chem. Phys. 1985, 82, 5732. (16) Carbajal-Tinoco, M. D.; Castro-Roma´n, F.; Arauz-Lara, J. L. Phys. Rev. E 1996, 53, 3745. (17) Larsen, A. E.; Grier, D. G. Nature 1997, 385, 230.

The concept of overcharging can be traced back to Kirkwood and Shumaker18 who showed that fluctuations in the net charge of a macroion lead to long-range attractions between the macroions. In their analysis they assumed that the charge fluctuations on each macroion were uncorrelated, which led in their series expansion of the pair potential to a first correction term proportional to 〈(∆Z1)2(∆Z2)2〉 ) 〈(∆Z1)2〉〈(∆Z2)2〉. They concluded that the magnitude of such charge fluctuations for electrically neutral species is comparable to that of permanent dipoles for certain biologically important molecules. It should be mentioned at the onset that there are other electrical mechanisms that may lead to attraction between macroions of like charge. These include short-range attractions due to anticorrelated charge density fluctuations between two rods, viz., 〈δF1δF2〉 < 0, as proposed by Oosawa19,20 with variations by Ha and Liu21 to explain the “bundling” of DNA, and by Rouzina and Bloomfield22 to include planar surfaces, depletion effects,23 and “bridging” by multivalent counterions.24-27 Also included in the long-range attraction category are the controversial Sogami-Ise28 effective pair potential with a long-range attractive tail, the “temporal aggregate” 29-32 and “temporal domain”33 models, and the “orbital model”34,35 based on the juxtaposition of potential field36,37 description of (18) Kirkwood, J. G.; Shumaker, J. B. Proc. Natl. Acad. Sci. U.S.A. 1952, 38, 863. (19) Oosawa, F. Biopolymers 1970, 9, 677. (20) Oosawa, F. Polyelectrolytes; Marcel Dekker: New York, 1971. (21) Ha, B. J.; Liu, A. J. Phys. Rev. Lett. 1998, 81, 1011. (22) Rouzina, I.; Bloomfield, V. A. J. Phys. Chem. 1996, 100, 9977. (23) Bhuiyan, L. B.; Vlachy, V.; Outhwaite, C. W. Int. Rev. Phys. Chem. 2002, 21, 1. (24) Linse, P.; Lobaskin, V. Phys. Rev. Lett. 1999, 83, 4208. (25) Linse, P.; Lobaskin, V. J. Chem. Phys. 2000, 112, 3917. (26) Linse, P. J. Chem. Phys. 2000, 113, 4359. (27) Hribar, B.; Vlachy, V. Biophys. J. 2000, 78, 694. (28) Sogami, I.; Ise, N. J. Chem. Phys. 1984, 81, 6320. (29) Schmitz, K. S.; Parthasarathy, N. In Scattering techniques Applied to Supramolecular and Nonequilibrium Systems; Chen, S.-H., Chu, B., Nossal, R., Eds.; Plenum: New York, 1981; p 83. (30) Schmitz, K. S.; Parthasarathy, N.; Vottler, E. Chem. Phys. 1982, 66, 187. (31) Schmitz, K. S. Biopolymers 1982, 21, 1383. (32) Schmitz, K. S.; Lu, M.; Gauntt, J. J. Chem. Phys. 1983, 78, 5059. (33) Fo¨rster, S.; Schmidt, M.; Antonietti, M. Polymer 1990, 31, 781.

10.1021/la0301590 CCC: $25.00 © 2003 American Chemical Society Published on Web 10/15/2003

Overcharging in Macroions

the distribution of counterions based on the potential field profile set up by the collection of macroions. Of particular relevance to the present study are the three-dimensional density profiles of both the counterions and co-ions in the presence of macroion clusters as determined by Brownian dynamics (BD) simulations.35,38 First, these profiles clearly indicate that the co-ions are expelled from the interior of macroion clusters whereas the counterions are drawn into these interiors. If the mechanism for overcharging pertains to macroions in a cluster, then these simulations indicate one need not consider co-ions in the surface simulations. Second, the “effective charge” as estimated from a “thermal radius”39-41 was found to vary with position of the macroion in the lattice, being virtually neutral at the lattice center.38 Third, the three-dimensional plots indicate that the interior counterions are not equally associated with the interior sides of participating macroions in the cluster.42 Unfortunately it could not be determined if those macroions with the higher local concentration of counterions could be classified as “charge reversal” macroions. However, the “near neutral” value of the “effective charge” of interior macroions of a cluster means that charge fluctuations will have an important contribution to the stability of the cluster. The present study expands upon the molecular dynamics (MD) studies of Messina and co-workers.10-13,43 in which the distribution of monovalent and divalent counterions for overcharging spherical macroions were considered. Of particular focus is the conclusion of Messina et al.13 that “... an isolated spherical colloid is always overcharged” and that the total electrostatic energy of a macroion can decrease to a certain minimum following addition of counterions in excess of what is required to neutralize it. It is worthwhile, therefore, to review briefly the salient features of their works relevant for our purposes. 2. Messina-Holm-Kremer Methodology on Macroion Charge Reversal One of the motivations of the Messina et al. studies10-13,43 was to investigate possible attraction mechanism between two like-charged macroions via metastable ionization. These authors began by showing how basic electrostatics and the Gillespie rule44,45 of physical chemistry may be combined to give the ground state of an isolated colloid (hard sphere with a point central charge) in the presence of counterions (electrons) on its surface: the ground state is overcharged with the counterions forming a regular geometric pattern maximizing their mutual separations on the sphere surface. The Messina and co-workers10,13,43 study is comprised of three parts: (1) molecular dynamics (MD) simulations for discretely placed macroion charges (denoted by the letters DCC); (2) MD simulations in which the macroion charge is “smeared” over the surface (denoted by the letters CC); and (3) an approximate theory based on the Wigner crystal (WC) framework, which yields an analytical (34) Schmitz, K. S. Handbook of Polyelectrolytes and Their Applications; Tripathy, S. K., Kumar, J., Nalwa, H. S., Eds.; American Scientific Publishers: Stevenson Ranch, CA, 2002; 3, 196. (35) Schmitz, K. S. Phys. Rev. E, in press. (36) Schmitz, K. S. Langmuir 1999, 15, 4093. (37) Schmitz, K. S. Phys. Chem. Chem. Phys. 1999, 1, 2109. (38) Schmitz, K. S. Langmuir 2001, 17, 8028. (39) Schmitz, K. S. Langmuir 2000, 16, 2115. (40) Sanghiran, V.; Schmitz, K. S. Langmuir 2000, 20, 7566. (41) Mukherjee, A. K.; Schmitz, K. S.; Bhuiyan, L. B. Langmuir 2002, 18, 4210. (42) Schmitz, K. S. Phys. Rev. E 2002, 66, 061403-061401. (43) Messina, R. Physica A 2002, 308, 59. (44) Gillespie, R. J. J. Chem. Educ. 1963, 40, 295. (45) Gillespie, R. J. Struct. Chem. 1998, 9, 73.

Langmuir, Vol. 19, No. 23, 2003 9601

expression for the purpose of characterization and employs a circular WC cell for the counterion-macroion interaction, and an approximation for the pairwise interaction of the excess (charge reversing) counterions. We review only the CC simulation method and the WC model. It should be mentioned that the DCC and CC simulations gave similar results. 2.1. CC Simulations. By use of Gauss’ theorem, the surface spread macroion charge was treated as a total charge at the center of the sphere. The counterions of valency Zc were assumed to interact with the charged sphere by a Coulomb potential and a Lennard-Jones potential to account for the excluded volume effects. Being a “soft core” potential, the latter allowed partial penetration of the counterions into the sphere’s surface. In applying the MD techniques at absolute zero temperature, that is, in strong Coulomb coupling, an important feature was the systematic annealing of the system after every fixed number of cycles to avoid traps of metastable states. Interesting findings were the facts that the counterions were seen to be highly ordered on the surface of the macroion and that the energy per counterion was inversely proportional to the average separation distance between the counterions on the macroion surface. This has relevance for the analytical model discussed below based on the WC ideas. 2.2. Analytical Model. The counterions were partitioned into two groups: the Nc neutralizing counterions and the n “overcharging” counterions. The reference state was taken to be the neutral particle composed of the macroion and the Nc neutralizing counterions. A “background charge” is defined in terms of a WC cell with a diameter equal to the average separation distance between the surface counterion charges and the surface concentration of counterions. Thus the average separation distance between the “surface bound” counterions is (A/ 4π(Nc + n))1/2 ) (1/(Nc + n))1/2R, where R is the macroion radius. The interaction energy of the N ) Nc + n counterions was referred to as the “correlation energy”. In addition to this interaction energy there is another term that can be taken to be an average pairwise interaction energy of the n overcharging counterions, where the number of such interactions was taken to be n2/2. Assuming that the interaction is Coulombic in origin, an approximation to the average separation distance between the n ions was assumed, for convenience, to also be the macroion radius R. Thus we define the energy gain ∆En relative to the reference state for (Nc + n) bound counterions in terms of their energy expression rewritten in our notation10,11

∆En ) -

2 RλeZc2 n2 Zc λe 3/2 3/2 (1) [(N + n) N ] + c c 2 R A1/2

where λe ) |e|2/4πor defines the charge parameter, r and o are the relative and vacuum permittivity, respectively, and R is a geometric dimensionless parameter related to the “background charge”. In the case of a planar surface, R ) 1.96. Obviously, a little error in R is introduced in the case of curved surfaces, because the area covered by the total number of WC cells is not exactly the same as the surface area of the spherical macroion. Note also that due to their Euler characteristic, the sphere does not possess a regular triangulation, and it can be shown that apart from a triangular lattice there have to be 12 defects with only 5 nearest neighbors. A good discussion of this and the relation of the minimization problem to the Thomson

9602

Langmuir, Vol. 19, No. 23, 2003

Mukherjee et al.

Table 1. Physical Parameters Used in the Present Calculations standard state temperature relative permittivity valence of counterion valencies of macroion ion counterion diameter radius of macroion

T ) 298 K r ) 16.0 Zc ) 2 Zm ) 50, 90, 180 σ ) 1.8 Å R ) 28.56 Å

problem have been given.46-48 Thus the counterions cannot be at exactly equidistant positions on the surface of the macroion. The more the ions try to be “equidistant” from each other (for a macroion of given radius), the closer the value of R is to the standard value. Thus R is a good indicator of some “uniformity” in ion distribution. The parameter R can be written in terms of ∆E1

R)-

((Zc2λe/2R) - ∆E1)A1/2 Zc2λe[(Nc + 1)3/2 - Nc3/2]

(2)

to obtain the operational expression

(

∆En ) -

)

Zc2λe [(Nc + n)3/2 - Nc3/2] n2 Zc2λe - ∆E1 + 2R 2 R [(Nc + 1)3/2 - Nc3/2] (3)

We shall refer to the above expressions (eqs 1-3) as the “WC equations” or the “WC theory” hereafter to distinguish the results obtained from the approximate analytical model and those from the simulations. 3. Models and Simulation Methods In this paper, we expand upon the works of Messina et al.13 for divalent counterions and explore an alternate simple technique based on electrostatic energy minimization (EM) and random sampling. To directly compare our EM methods with the MD results of Messina et al.,13 we likewise employed a sphere with radius R ) 28.56 Å and limit the counterion charge to the value Zc ) 2 with diameter σ ) 1.8 Å for all the calculations. The model systems in our simulations are comprised of an isolated macroion with bare charge -Zm|e| surrounded by N ) Nc + n number of small counterions with charge q ) Zc|e|, where Nc is the number of neutralizing counterions, viz., Zm|e| ) -NcZc|e| is the neutral state taken to be the reference state, and n is the number of overcharging counterions. The parameters used in our study are summarized in Table 1. 3.1. Spherocylinder and Spheroid Shapes. We extend these studies to the nonspherical shapes of spherocylinders and spheroids. The spherocylinder is represented by two hemispheres connected by a cylinder of the same radius as that for the hemispheres, whereas a spheroid may be thought of as arising from a continuous distortion of a sphere by either compressing (oblate spheroid) or stretching (prolate spheroid) the sphere from two opposite points. Thus the choice of geometries represents topologically equivalent shapes. Another reason to study these geometric shapes is that many biological organisms and chemical macromolecules have similar shapes, viz., globular proteins. The different geometries are identified by a “shape factor” τ (see Appendix A), which is defined as τ ) L/r, where L is a characteristic length that measures the (46) Thomson, J. J. Philos. Mag. 1904, 7, 237. (47) Perez-Garrido, A.; Moore, M. Phys. Rev. B 1999, 60, 15628. (48) Bowick, M.; cacciuto, A.; Nelson, D. R.; Travesset, A. Phys. Rev. Lett. 2000, 89, 185502.

Figure 1. Spherocylinder. A and B are the centers of the two hemispherical end caps which are L distance apart (L is also the length of the cylinder).

deviation from spherical shape and r is a characteristic radius. A sphere is therefore defined by L ) 0 and the parameter r is identified with the radius of the sphere. The construction of the spherocylinder may be thought of as cutting a sphere into two equal halves and then separating their centers by a distance L. The values of L, and therefore τ, are constrained to be positive. If the spheroid is obtained by compressing the sphere, then the values of L, and therefore τ, are negative. The resulting structure is the oblate spheroid, or the shape of a discus. If the spheroid is obtained by stretching the sphere, then the values of L, and also τ, are positive. The resulting spheroid is a prolate spheroid. Since the spherocylinder and prolate spheroid are similar, calculations are not performed in the present study on prolate spheroids. The geometry and parameters are shown in Figures 1 and 2. The length of the cylinder is L ) rτ. In the spheroid, the centers A and B of the curvatures of the two hemispherical sections A and B are at a distance L apart. Obviously, if τ ) 0, the two nonspherical geometries reduce to sphere. Rather than perform our simulations with a fixed surface charge density, we have chosen to fix the volume charge density. The volume charge density as calculated for the sphere F0 ) 3Zm/4πR3 is fixed for all geometries considered herein. This choice is made to ensure that the total macroion charge, -Zm|e|, of the topological set of shapes remains constant. For example in the case of a spherocylinder, when the length L is changed, the radius of the hemispherical end caps r is also changed to keep same volume. Although τ can assume any value, we have considered -1 e τ e +2 in this paper. We consider the following assumptions for the internal charge distributions of nonspherical macroions: in a spherocylinder, a line charge of length L along the axis of the cylinder carries the total charge -Zm|e| as shown in Figure 1; and in a spheroid, the total charge -Zm|e| is assumed to be a point charge and is situated at the midpoint C of the line joining the centers A and B as shown in Figure 2.


Langmuir, Vol. 19, No. 23, 2003 9603

interaction of the surface located counterion with the spheroid charge is

(

Figure 2. Spheroid. A and B are the centers of the curvatures of the two equal sections a and b of a sphere of radius r. Centers A and B are at L distance apart. The point O is considered as the center of the spheroid where the total charge Zm of the macroion is assumed to be concentrated.

3.2. Energy Considerations. The simulations take into consideration the pairwise interaction between all of the N ) Nc + n counterions. Thus within the framework of the WC theory the energy difference in the Messina et al.13 notation and in the notation we have is ∆En ) EN ENc where Nc ) Zm/Zc. It is important to note that the dielectric constant of the macroion is always assumed to be the same as that of counterions and the supporting medium. This is a rather common assumption made in the literature to avoid the problem related with electrostatic image due to the discontinuity of dielectric constant at the interface. It ought to be mentioned however, that effects of different dielectric constants and charge discretization have been treated by Messina in some very recent work.43 The Coulombic potential energies EN for the three geometries can be written straightforwardly as the general expression for the counterion on the surface

(∑ ) Zc2

EN ) -Ecs + λe

i 0): Line Charge. The cylindrical section of the spherocylinder has a radius r and a length L. The expression of the counterion interaction with a line charge on a spherocylinder is quite complex. The expression depends on whether the counterions lie either on the hemisphere or on the cylindrical part of the spherocylinder. If the z-coordinate of the counterion is zi > L/2, then it lies on the upper hemisphere, whereas if zi < -L/2, then it lies on the lower hemisphere. The general expression for the energy of interaction of the counterions with the spherocylindrical macroion is

Ecs ) λe

ZmZc L

N

(

)

{1 - sin(ψ2i)} cos(ψ1i)

loge ∑ {1 + sin(ψ i)1

1i)}

cos(ψ2i)

(7a)

where

ψ1i ) tan-1

(

)

(((L/2) - zi)

ψ2i ) tan-1

r cos(θi)

(

)

(L/2) + zi r cos(θi)

(7b)

If the counterion is on the cylindrical part of the spherocylinder, then θi ) 0 and the positive sign in the first of the eq 7b applies. If the counterion is on either of the hemispherical endcaps, the negative sign is to be taken and

θi ) sin-1

(

)

zi - (L/2) r

(7c)

In the particular case when zi ) (((L/2) + r), it is convenient for numerical purposes to utilize the expression

ZmZc r loge Ecs ) -λe L r+L

(

)

(7d)

(iv) Spherocylinder (τ > 0): Point Charge. In this case the interaction energy Ecs is calculated in the same way as for a spheroid, as given by eq 6 where the charge Zm is located at the origin of the cylinder coordinate system. 3.3. Simulations. Since ground-state energy configuration is defined at absolute zero temperature,13 there is no random motion of the associated counterions. Therefore it is not possible to apply the standard procedures of energy selection of the Metropolis algorithm of Monte Carlo simulations in this study. Moreover, only one final configuration is required for this system which has the lowest possible total electrostatic potential energy.

9604


In this technique, to achieve that final configuration we followed the following two stages. We first established the configuration of the neutral state where the Nc counterions exactly neutralize the macroion. The counterions were distributed on the surface of the macroion by giving them random coordinates. An ion was then chosen at random and was moved randomly on the macroion surface. The total electrostatic energy of the whole system was then calculated. The move was accepted only if it caused a decrease in total energy. In all cases, each ion attempted approximately one million moves before the final configuration of the system was reached. The second stage is that for overcharging the macroion. Starting from the neutral state, one counterion was added to the system and the energy was then minimized following the above procedure. The displacement parameter, which is defined as the extent of maximum movement of an ion, adjusted itself continuously (in a random manner) until an ion move was accepted. This actually helped to avoid the trap of metastable states (local minima) and thus ensured the maximization of mutual ion distances very quickly leading to minimization of the interaction energy. Then another counterion was added and the procedure repeated. This procedure was followed for each addition. Specifically, the continuous adjustment of the displacement parameter was done in the following manner. At the very beginning of the simulation run, a displacement parameter was assigned which was approximately equal to the average distance between two adjacent ions. By use of this displacement parameter the system was brought to its first metastable state and no further decrease in energy was found. In such a situation, the displacement parameter was changed when the energy was seen to decrease again until the system reached the next metastable state. Actually, when the system has reached a metastable state, no movement of ions is accepted since this does not decrease the energy. To break this special arrangement of ions, a different extent of ion displacements is needed. Normally, during the EM run, the displacement parameter gradually decreases while the system passes through a series of metastable states of increasingly lower energy. In general for all systems considered in this work, the displacement parameter reduced to less than 10-6 Å within one million moves per counterion and no more decrease in energy occurred in some particular cases. But in most of the cases we have found (after a sufficient number of iterations) that even after the displacement parameter has adjusted to a value less than 10-6 Å the energy continues to decrease as the simulation proceeds, but the energy changes are infinitesimally small. In such cases we assume that the system is very close to the true ground state. We remark further that for the purposes of the present study, knowledge of the exact ground state energy is really not necessary since one wants only to be sufficiently close to the ground-state energy. The computer routine has been made in such a way that it can calculate minimum energy configurations for any of the three models used in this study by choosing only one of the three following values of the “shape factor”: τ ) 0 (sphere), τ < 0 (spheroid), τ > 0 (spherocylinder), respectively. 4. Simulation Results As indicated earlier, keeping in conformity with the previous work of Messina et al.13 we have considered the same parameters associated as used by these authors for ease of comparison. We also represent the energies as

Mukherjee et al.

Figure 3. Ground-state electrostatic energy of a charged spherical macroion as a function of overcharging counterions. Open circles are from Messina,10,13 filled circles are from this work, and the solid lines are from eq 3. The neutral state energies of a spherical macroion of corresponding charges have been taken as the origin of the plot. Curves a, b, and c are for Zm ) 50, 90, and 180, respectively.

being “normalized” to the thermal energy under standard thermodynamic conditions, viz., T ) 298 K. The standard state thermal energy is kBT ) 4.1124 × 10-21 J/particle. 4.1. Results for Spherical Macroion. In Figure 3, we have plotted the difference between the ground-state configuration energy and the reference energy of a spherical macroion versus the number of overcharging counterions. The solid circles are for the present study and for comparison we have included the results of Messina et al.13 shown by the open circles. These energies are “normalized” to the thermal energy at the standard state temperature, kBT ) 4.1124 × 10-21 J/particle. These data show a clear minimum in the energy profile. We denote the reduced energy minimum as ∆Emin/kBT, which occurs at nmin ()n(∆E/kBT ) ∆Emin/kBT)), and we list some representative values in Table 2. For the lowest macroion charge corresponding to Zm ) 50, the energy variation for the present study is almost the same as those given in Figure 5 of Messina and coworkers.13 With the increase of macroion charge, the results clearly show a gradual deviation from their corresponding results. The reason for this discrepancy is probably due to the choices of the potentials in the two studies. Messina and co-workers employed the LennardJones (LJ) potential as a “soft potential” for the excluded volume interactions whereas the hard-core potential was used in the present work. The qualitative difference in these two studies may be explained as follows. The LJ potential allows the counterions to penetrate up to a distance 0.5σ into the macroion surface. One would therefore expect the Messina curve to be lower than that of the present work. However, the repulsive counterioncounterion interaction would likewise be stronger in the Messina case since their interparticle distance would also be smaller than in our simulations. Thus the larger number of repulsive pairwise interactions overrides the anticipated increase in the counterion-macroion interaction and the Messina curves thus lie above ours as the number of overcharging counterions is increased. In case of smaller systems (Nc ) 25), however, the counterions


Langmuir, Vol. 19, No. 23, 2003 9605

Table 2. A Comparison between the Ground State Electrostatic Energies Zm

n

EMa

WCb R ) 1.956

WCc R ) 1.92

WCd R ) 1.96

MDe

50

0 1 2 3 4 5 6

0 -18.060 -31.955 -40.843 -45.051 -45.192 -40.651

0 -18.060 -31.611 -40.644 -45.152 -45.129 -40.569

0 -17.683 -30.848 -39.489 -43.598 -43.169 -38.195

0 -18.102 -31.695 -40.772 -45.324 -45.346 -40.832

0 -17.988 -31.747 -40.423 -44.348 -44.159 -39.206

Zm

n

EMa

WCb R ) 1.935

WCc R ) 1.88

WCd R ) 1.96

MDe

90

0 1 2 3 4 5 6 7 8

0 -24.596 -44.747 -61.233 -71.814 -78.566 -80.346 -77.389 -69.998

0 -24.596 -44.578 -54.945 -70.693 -76.814 -78.320 -75.193 -67.436

0 -23.884 -43.148 -57.788 -67.801 -73.185 -73.937 -70.053 -61.532

0 -25.005 -45.402 -61.187 -72.357 -78.910 -80.843 -78.152 -70.834

0 -24.379 -44.220 -60.237 -70.302 -76.464 -77.567 -73.859 -65.639

WCb

WCc

WCd

Zm

n

EMa

R ) 1.93

R ) 1.88

R ) 1.96

MDe

180

0 1 2 3 4 5 6 7 8 9 10

0 -36.132 -67.980 -93.616 -116.330 -133.502 -146.658 -154.154 -157.432 -156.793 -149.601

0 -36.132 -67.569 -94.308 -116.349 -133.699 -146.333 -154.271 -157.511 -156.045 -149.875

0 -34.692 -64.679 -89.962 -110.540 -126.409 -137.571 -144.024 -145.766 -142.798 -135.117

0 -36.272 -67.850 -94.731 -116.914 -134.399 -147.185 -155.270 -158.653 -157.333 -151.310

0 -35.270 -66.340 -92.110 -113.160 -129.480 -141.040 -147.580 -149.560 -146.380 -138.710

a Simulations, this study. b From eq 1 using the indicated value of R, this study. c From Messina et al.13 using eq 1 and the indicated value of R. d From eq 1 using the standard value of R, this study. e Molecular dynamics results of Messina et al.13

seldom encounter each other or penetrate into the macroion due to comparatively weaker Coulombic interaction. Thus results of our study for smaller systems are very close to those from MD calculations.13 In larger systems (Nc ) 45, 90) however, the LJ potential in MD calculation plays an important role in ion distribution due to shorter distances (e21/6σ) between the counterions and to their being able to penetrate inside the macroion. Obviously, for a higher number of counterions, the two techniques (MD and EM) actually consider different physical conditions and are therefore not directly comparable. We also remark in passing that the result out of the analytical model, namely, eqs 1-3, shown by the solid line in Figure 3 is in close correspondence with the simulations. In the application eq 3, the quantity R was calculated from eq 2 using the simulation result for the first overcharge ∆E1. 4.2. Results for Nonspherical Macroions. The excess energy for different nonspherical shapes as a function of overcharging counterions is shown in Figure 4 for Zm ) 50. We have studied other macroion charges and found similar results. The lowest curve is for sphere (which is the same as the first curve of Figure 3). Each curve is labeled by corresponding factor τ with τ ) 0, τ > 0, and τ < 0 indicating sphere, spherocylinder, and spheroid, respectively. In both the nonspherical cases, the minima of the interaction energy are increasing (becoming more positive) with the increase of |τ|. The lowest value of minima among all geometries occurs for the case of spheroidal macroion with 0 > τ g -0.1 (this curve lies very close to the curve of the sphere and is therefore not

Figure 4. Ground-state electrostatic energy of charged macroions of three different geometries and different τ as a function of overcharging counterions n. Zm ) 50 and the geometries are the sphere (τ ) 0), spheroid (τ < 0), and spherocylinder (τ > 0). The neutral-state energy of each size and geometry is the origin of corresponding curves. For both the nonspherical geometries the lowest energy points are shifting upward with the increase of |τ|. The dotted lines are polynomial (third degree) fits (except the line for τ ) 0, which is from eq 3) to guide the eye.

shown in Figure 4). We have conducted a systematic study of spheroids and spherocylinders with τ ranging from -1 to 2 and found that in no case does the minimum energy become lower than, or even become equal to, that for a spheroid with τ ≈ -0.1. In both cases these curves bend upward (toward positive) with the increase of |τ|. From Figure 1 and Figure 2, it is obvious that increase of |τ| for a spheroid means decrease of its thickness, while for spherocylinder an increase of |τ| implies an increase of its length. Since in both geometries the volumes are kept always the same as that of a sphere of radius R ()28.56 Å), the increase of |τ| makes them narrower, bringing the surfaces closer to the line charge for spherocylinder and to the central point charge for a spheroid. Due to this, the interaction energy should decrease (becomes more negative) but we see the opposite in Figure 4 and in Figure 5. This can be explained in the following way: In a spheroid, with the increase of |τ| the counterions come closer to one another causing increased repulsion among them. The repulsive part of the total interaction energy becomes more significant than the attractive one. In a spherocylinder, increase in length decreases the line charge density as 3Zmr2/4πR3 and consequently the attractive interaction energy decreases as 3ZmZcr2/16π2orR3. But in the spherocylinder case an increase of surface area also decreases the repulsion among the counterions. Due to this competitive decrease of attractive and repulsive interaction energies, increase in total energy with |τ| is much slower for a spherocylinder than that for a sphere. In Figure 5 we have plotted maximum energy gain due to overcharging against τ. The maximum energy gain is the energy difference between the neutral state and the maximally overcharged state (i.e., depth of the minimum) for any τ. This figure also clearly depicts that the energy gain is the lowest for spheroids of 0 > τ g -0.1. It is seen that for a system with divalent counterions numbering around 210 (Zm ) 420, τ ) -0.01) and around 175 (Zm ) 350, τ ) -0.10), the energy gain for a spheroid is bigger than that of a sphere, while for higher numbers of counterions the roles reverse. The energy gain is increasing

9606


Mukherjee et al.

Figure 5. The maximum energy gain due to overcharging for microions of all the three geometries and different charges as a function of τ. The fit is a third degree polynomial for τ < 0 and line for τ > 0 to guide the eye.

Figure 7. Minimum energy configurations of (A) 60 and (B) 92 divalent counterions on a spherical macroion (neutral). The counterion positions are shown by dots and circles. The dots and solid lines are on the front and the circles and dashed lines are on the backside of the macroion sphere. Figure 6. Neutral-state energy of macroions with different geometries and bare charges as a function of τ. The neutralstate energy of the sphere is taken as origin. The dotted lines are third-order polynomial fits to guide the eye.

linearly with the increase of τ for a spherocylindrical macroion, while for a spherical macroion it is nonlinear. At lower Zm and τ it fluctuates around some fitted line. With the decrease of τ, the surface of the spheroid comes closer to the point C in Figure 2 as stated above and this causes a very strong attractive interaction between counterions and macroion. By the same token the area where the counterions find themselves shrinks and repulsion among them also becomes higher. In such a situation any excess counterion, due to overcharging, would further enhance both attraction and repulsion. Such competing effects are likely to cause fluctuation in total energy. This is clearly depicted by the curve labeled τ ) -0.75 of Figure 4 where addition of the fourth excess ion increases the energy than from that for the third excess ion, and further addition (the fifth) decreases the energy again. The expected energy gain (the fitted line) is not maintained.

Figure 6 shows that for increase of |τ| the neutral state energy of both geometries decreases but the decrease is much more rapid for a spheroid than that for a spherocylinder. The most significant feature of Figure 6 is that the neutral state energy of a spherocylinder also has a minima at some length of the cylinder similar to the overcharging curves in Figure 3. It is obvious from Figure 6 that the energy of neutral state is a continuous function of τ which shows the dependence of Coulombic interaction energy with shape and size of a charged particle. In nonneutral cases we also found curves similar to that in Figure 6. 4.3. Geometry and Counterion Distribution. Shown in Figure 7 are divalent counterion distributions for neutral macroions of charges Zm ) 120 and 180. The counterions distribute themselves as far apart as possible to minimize the energy, giving rise to pentagonal and hexagonal patterns with an underlying triangular geometry. The counterion distributions for the spherocylinder and spheroid are shown in Figure 8. In contrast to the results for the sphere, the counterion distributions for these two geometries are not uniform on the surfaces.


Langmuir, Vol. 19, No. 23, 2003 9607

Figure 8. Minimum energy configurations. The values for τ are (0.25, (0.5 and (0.75. In the case of a spheroid (τ < 0), the distributions are shrinking with increasing |τ| and are gradually losing regular well-defined arrangement. The distribution pattern of a spherocylindrical macroion (τ > 0) for lower τ ()0.25) is very similar to that of a spherical one (Figure 7) and gradually assumes another well-defined pattern with three distinct lines at the middle of the cylinder with increasing |τ|. The notation of the symbols is the same as that given in Figure 7.

9608


To explain the counterion distributions in Figures 7 and 8, one must be reminded that they are based solely on the minimization of the electrostatic energy that represents the ground state of these complexes. The pattern of the counterion distribution on the surface of a macroion is therefore a direct consequence of both the geometry of the particle and the macroion charge distribution. Of course the repulsive interaction between the counterions, both the excluded volume and electrostatic contributions, plays a role as well. What is important is the relative contributions to the interaction energy from the counterion-macroion and the counterion-counterion interactions. Consider first the spherical geometry. Because of the placement of the macroion charge at the center of the sphere, all of the counterions will experience the same interaction with the macroion regardless of their positions on the surface of the sphere. Hence the distribution of the counterions on the sphere is due to the repulsive interactions between the counterions themselves. The sphere with a central charge has no influence on the counterion distribution. Now consider the shape with the macroion charge represented as a point charge at the origin of the descriptive coordinate system. Counterions at the equator of the spheroid will experience a weaker interaction with the macroion than those counterions at the polar caps of this geometric shape. Hence the “flatter” the sphere, the more likely it is to find counterions in the polar regions with an absence of counterions at the equator. If the counterions were point charges that did not interact with each other, but only with the macroion, then one would expect all of the counterions to be located at either pole location. This is a direct result of energy minimization. However, the counterions do interact with each other and this repulsive interaction thereby pushes the counterion distribution toward the equator. Had we chosen to represent the internal charge distribution of the macroion as a ring in the interfacial plane of the two halves, then energy minimization would direct the counterions to congregate at the location of this charge ring and the mutually repulsive interactions would then spread the distribution toward the poles as well as the equator. We now address the charge distribution within the spherocylinder. We represented the macroion charge as a line charge in our simulations, which is the usual representation for cylindrical particles. The Manning condensation model for DNA49 is one such example. Following the above arguments, one would expect the counterion distribution to be more or less uniform on the cylindrical surface with a less dense distribution in the cap regions. Had we represented the macroion charge as a point at the origin of the descriptive coordinate system, then we would expect the counterion distribution to be maximum at the waist of the cylinder and perhaps no counterions to be found at the cap regions. This is one of the reasons for representing the spherocylinder charge as a line charge rather than a point charge, for to our knowledge counterions are distributed throughout the surface of cylindrical macroions such as DNA. These arguments are borne out in Figures 9 and 10. For example, in Figure 9 the deviation in excess energy in a spherocylinder for a point charge relative to a line charge increases with τ as expected, while in Figure 10 (for a point charge) the dense distribution of the counterions around the central region of the spherocylinder relative to that in the endcaps may be noted. (49) Manning, G. S. Acc. Chem. Res. 1979, 12, 443.

Mukherjee et al.

Figure 9. Comparison of point versus line charge representations for spherocylinders. Filled symbols are for line charge, while the open symbols are for point charge. All the results correspond to Zm ) 50. The dotted lines are third-order polynomial fits to guide the eye.

Figure 10. Minimum energy configuration of 90 divalent counterions on a spherocylinder for τ ) 2 and a point (macroion) charge (Zm ) 180) at the origin of the cylindrical coordinate system. The notation of the symbols is the same as that given in Figure 7.

5. Discussion One of the reasons for studying the properties of overcharged systems is to provide an explanation for the heterogeneous two-state structure of colloidal systems and the “bundling” of rodlike particles such as DNA. Indeed Messina, Holm, and Kremer11 have shown that for divalent counterions two unlike charged spheres do show an attraction at short distances. Just as in the case of mass accretion in a binary star system, the counterions from the lesser charge sphere are drawn to the surface of the more highly charged sphere to form some kind of “ionic bond”, to use the terminology of chemical bond theory. This finding supports the study of Barouch and Matijevic´50 in which spheres of different size exhibit short-range attraction due to induced surface charge while at long distances show a repulsive interaction. The problem, however, is that monovalent and not divalent counterions are present for the experimental systems in which heterogeneous structures are observed, and then under low ionic strength conditions. One must therefore ask, under what conditions might a “charge reversal” mech(50) Barouch, E.; Matijevic´, E. J. Chem. Soc., Faraday Trans. 1985, 81, 1797.


Langmuir, Vol. 19, No. 23, 2003 9609

Table 3. Comparison of the WC Parameter r Obtained from EM and MD Methods

Table 4. Neutral State Energy and Minimum Energy for a Spherical Macroion at Selected Values of Zm

Zm

Nc

Ra

Rb

Zm

EZma/kBT

Emina/kBT

nmin

4 6 8 10 20 30 32 50 90 128 150 180 288 360

2 3 4 5 10 15 16 25 45 64 75 90 144 180

1.894 1.884 1.894 2.005 1.911 1.954 1.956 1.956 1.935 1.956 1.990 1.930 1.867 1.835

1.890 1.970 1.920 2.020 1.930 1.910 xxxx 1.920 1.880 xxxx 1.910 1.880 xxxx 1.860

50 90 180

-1871.027 -5786.107 -22189.829

-1916.219 -5866.453 -22347.261

5 6 8

From eq 2 using simulation results, this study. b From Table 2 of Messina et al.13 using MD results.

a

In units of kBT ) 4.1124 × 10-21 J/particle.

Table 5. Neutral State Energy and Minimum Energy for a Nonspherical Macroion at selected values of Zm and τ τ 0.25 +0.50

a

anism be a candidate, in general, for the explanation of grouping of macroions in real systems. Recent Brownian dynamics (BD) simulations of macroion clusters may provide an answer: Using the concept of a “thermal radius” to define an “effective charge” of the macroion,40,41 it was determined that the “effective charge” of a macroion at the center of a “diamond-shaped” cluster of seven spheres was about 8% of its “bare” charge.38 More recently BD simulations showed that the co-ions were virtually excluded from the interior of a macroion cluster while the counterions were drawn into the interior of the cluster.35 5.1. Comparison of the EM and MD Results. Our comparison with the WC model lies solely with the parameter R as calculated from the first overcharging counterion in accordance with eq 2. When the macroions of nonspherical shapes are gradually overcharged in the same way as done in Figure 4, the potential energy is expected to vary similarly with excess counterion number irrespective to the shape and size of the macroions. In Table 2 we summarize the numerical results for direct comparison with the corresponding numerical values from the Messina et al. work13 (Messina, personal communication). Columns 3, 4, and 7 have already been shown graphically in Figure 3, where the dots are column 3, solid lines are column 4, and the open circles are column 7. Using first overcharge ∆E1 from column 3, WC parameter R has been calculated from eq 2. Applying those R values in eq 1, we have calculated column 4. The close proximity between the two columns, 3 and 4, clearly depicts excellent agreement between this technique and the WC theory. We have calculated column 5 and 6 using R given by Messina et al.13 and WC standard R )1.96, respectively, from eq 1. Since in all the three cases (Zm ) 50, 90, 180) the values of R of our results are very close to the value of 1.96, the results are also very similar (cf. columns 3 and 6). Column 7 is the MD simulation result of Messina et al.13 All the results in Table 2 are in good numerical agreement overall testifying to the validity and accuracy of the present technique. In Table 3, we have calculated R for various cases and compared with those given in Messina et al.13 Again we see good agreements between the two sets of R calculated using two different techniques. Furthermore, we have seen that all data points are well fitted with eq 3 based on a simple version of the WC theory.13 The extension of this idea to other macroion geometries shows that the groundstate energy vs overcharge curves have similar shapes for all the three geometries, which indicates that ionic correlations among counterions over the surface of the macroions are similar to those in the spherical case. In

+2.0

-0.25 -0.5 a

Zm

EZma/kBT

Emina/kBT

nmin

50 90 180 50 90 180 50 90 180 50 90 180 50 90 180

-1879.833 -5799.875 -22222.134 -1895.727 -5847.746 -22400.567 -1945.764 -5999.898 -22955.570 -1895.171 -5849.472 -22426.462 -1946.168 -6254.188 -22946.247

-1922.658 -5878.506 -22376.882 -1937.122 -5922.952 -22548.667 -1981.000 -6053.559 -23076.190 -1938.312 -5930.102 -22580.130 -1986.935 -6323.534 -23091.701

4 6 8 4 6 8 4 6 8 4 6 8 4 5 7


Table 6. Neutral State Energy and Minimum Energy for a Spherocylindrical Macroion at Selected Values of Point and Line Charges Zm and τ τ line charge

+0.25 +0.5 +1.0 +2.0

point charge

+0.5 +1.0 +2.0

a

Zm

EZma/kBT

Emina/kBT

nmin

50 90 180 50 90 180 50 50 90 180 50 90 180 50 90 180 50 90 180

-1879.833 -5799.875 -22222.134 -1895.727 -5847.746 -22400.567 -1936.196 -1949.764 -5999.898 -22955.570 -1907.955 -5879.958 -22639.546 -2003.410 -6150.542 -23494.611 -2217.915 -6807.293 -25907.698

-1922.658 -5878.506 -22376.882 -1937.122 -5922.952 -22548.667 -1972.303 -1981.000 -6053.559 -23076.190 -1947.873 -5954.094 -22663.630 -2037.369 -6214.055 -23623.004 -2248.931 -6846.188 -25990.179

4 6 8 4 6 8 4 4 7 7 4 5 7 4 6 7 3 4 6


Tables 4 and 5 numerical values of the neutral state energy and minimum energy for a spherical and a nonspherical macroion, respectively, are presented at selected values of Zm and τ. The corresponding results for a spherocylinder with a point (macroion) charge are shown in Table 6. In this study we have not examined WC theory for nonspherical geometries. Still, in this case, the WC theory is expected to provide negligible errors in calculating the electrostatic potentials for small |τ|, since the surface area changes only slightly within this range of τ. For example, it can be shown that

{

Sx - Ssphere 0.086, τ ) 2 ) 0.087, τ ) -1 Ssphere

(8)

where Sx is the surface area of a nonspherical macroion. 5.2. Analytical Expression for EM Results. There is a clear difference between the mathematical description

9610


Mukherjee et al.

of the simulations presented in this study and the WC model. In the simulations the pair interactions between all of the counterions are taken into explicit consideration (cf. eqs 4-7). This means that there are N(N - 1)/2 counterion pair interactions in the simulation. In contrast the WC model explicitly considers the pair interactions between the “overcharging” counterions leading to the number of n(n - 1)/2 counterion pair interactions. It is noted that eq 1 has the term n2/2 instead of n(n - 1)/2. This is because the leading term absorbs the difference in R, i.e., the Zc2λe/2R term in eq 2. In our proposed model we count all of the pair interactions by a modification of the Scatchard approach51 which employs average interactions. In this case we write for the average energy 〈EN〉 of the macroion-counterion system as

(

〈EN〉 ) -ZmZc〈A〉N + Zc2〈B〉

)

N(N - 1) 2

) -ZmZc〈A〉(Nc + n) + Zc2〈B〉

(9)

(

)

(Nc + n)(Nc + n - 1) 2

where 〈A〉 and 〈B〉 represent some average function of distance for energies of interaction of the counterion with the macroion (the 〈A〉 term) and the counterions with other counterions (the 〈B〉 term). To compare with the WC expressions, we expand this expression to arrive at the quadratic expression in n

〈EN〉 ) g0 + g1n + g2n

2

(10)

where

g0 ) -ZmZc〈A〉Nc + Zc2〈B〉

(

)

Nc(Nc - 1) 2

(11)

corresponds to the reference energy (neutral macroion)

g1 ) Zc2〈B〉

(

)

(2Nc - 1) - ZmZc〈A〉 2

(12)

Zc2 [2Nc(〈B〉 - 〈A〉) - 〈B〉] ) 2 where we used the identity Nc ) Zm/Zc, and

g2 ) Zc2〈B〉/2

(13)

We can now interpret the shapes of Figures 3 and 4. Let us write the energy difference represented in these figures in terms of 〈EN〉

∆En ) 〈EN〉 - 〈ENc〉 ) g1n + g2n2

(14)

At low numbers of overcharging counterions an increase in the value of n results in a decrease in the value of Etotal only if g1 < 0 which means that 〈A〉 is greater than 〈B〉 (2Nc - 1)/(2Nc). It is anticipated that 〈A〉 is independent of the number of counterions since this term reflects the average interaction of the counterion with the macroion. However, the term 〈B〉 is dependent upon the number of counterions N. This is because the average pairwise interaction is dependent upon the average spacing between counterions. (51) Scatchard, G. Ann. N. Y. Acad. Sci. 1949, 660.

As the number of counterions on the surface increases, the average spacing between them decreases and thus 〈B〉 increases. Therefore we may write for the derivative of eq 14

( )

( ) ( )

∂(∆En) ∂g1 ∂g2 2 ) g1 + n + 2g2n + n ∂n ∂n ∂n

(15)

Zc2n ∂〈B〉 2 ∂N

) g1 + 2g2n + (2N - 1)

where the identities (∂gj/∂n) ) (∂gj/∂〈B〉)(∂〈B〉/∂N)(∂N/∂n) ) (∂gj/∂〈B〉)(∂〈B〉/∂N) and N ) Nc + n were used. If the number of surface counterions is large, one may make the approximation (∂〈B〉 /∂N) ) 0. The physical reason behind this approximation is that the inclusion of one more counterion on a surface with a large number of counterions will not significantly reduce the average distance between the counterions. This approximation is within the spirit of the Scatchard approach.51 With this further approximation the location of the minimum in the ∆En profile is

nmin ) -

g1 2g2

) Nc

[

(16)

]

[

]

〈A〉 1 Zm 〈A〉 1 -1 + ) -1 + 2 Zc 〈B〉 2 〈B〉

Since nmin must be a positive number if a minimum exists, it follows that either 〈A〉 is greater than 〈B〉 or else 0 g Nc[(〈A〉/〈B〉) - 1] > -(1/2). To further interpret the results of our simulations with the present analytical model, we express the location of ∆Emin/kBT to the total number of counterions Nmin ) Nc + nmin, viz.

Nmin )

1 Zm 〈A〉 + 2 Zc 〈B〉

(17)

We use the data in Tables 4-6 to estimate a value for the ratio 〈A〉/〈B〉 on the basis of eq 17. The results are summarized in Table 7. It is a curious feature of this analysis that the values of the ratio 〈A〉/〈B〉 are virtually independent of the shape and charge distribution of the macroion even though their actual distribution of the surface counterions differs greatly. It may thus be inferred that the relative distance parameters dictate the location of Nmin rather than the absolute distances. The physical reason for this observation is that the counterions are drawn to the most favorable region governed by the counterion-macroion interaction (smallest distance) thus also resulting in a smaller average separation distance between the surface adsorbed counterions. It is noted that 〈A〉 is dependent only on the characteristics of the macroion, viz., its shape and charge distribution. In contrast 〈B〉 is dependent on both the properties of the macroion and also the number of counterions. As the number of counterions increases, the value of 〈B〉 will increase because their separation distance decreases. In a crude way one would therefore expect 〈B〉 to be inversely proportional to the square root of the area delegated to each particle, i.e., 〈B〉 ≈ (N/A)1/2. In this regard the present analytical model, which is tailored to the energy expression use in our simulations, differs substantially from the WC theory. According to the WC expression of eq 1, the distance dependence of the n2 term is the radius of the macroion sphere and thus independent


Langmuir, Vol. 19, No. 23, 2003 9611

Table 7. Relative Average Energies 〈A〉/〈B〉 at Nmin macroion spheroid

sphere

spherocylinder line charge

spherocylinder point charge

a

τ

Zm

nmin

Nmin

〈A〉/〈B〉a

-0.5 -0.5 -0.5

50 90 180

4 5 7

29 50 97

1.14 1.10 1.07

-0.25 -0.25 -0.25

50 90 180

4 6 8

29 51 98

1.14 1.12 1.08

0.0 0.0 0.0

50 90 180

5 6 8

30 51 98

1.18 1.12 1.08

+0.25 +0.25 +0.25

50 90 180

4 6 8

29 51 98

1.14 1.12 1.08

+0.5 +0.5 +0.5

50 90 180

4 6 8

29 51 98

1.14 1.12 1.08

+2.0 +2.0 +2.0

50 90 180

4 7 7

29 51 97

1.14 1.14 1.07

+0.5 +0.5 +0.5

50 90 180

4 5 8

29 50 98

1.14 1.10 1.08

+1.0 +1.0 +1.0

50 90 180

4 6 7

29 51 97

1.14 1.12 1.07

Calculated from eq 17.

of the number of counterions on the surface of the sphere. A detailed comparison of the WC and our modified Scatchard model will be presented elsewhere. 6. Concluding Remarks In this paper, we have presented a simple technique based on electrostatic Coulombic interaction energy minimization and seen its performance vis-a`-vis a published MD study13 in some special cases. Application of the technique to different geometries suggests it to be efficient, straightforward to use, and economic in terms of computer time usage. Note that, in the case of MD simulation, the system needs to be heated and cooled periodically13 to avoid the local minima, while in the present technique this is done automatically through selfadjustment of the displacement parameter. This is not to say that the technique will work under all physical circumstances. For instance, it is likely to be difficult in situations where large global rearrangements of the counterions are accompanied by a change of step length or the displacement parameter. However, for systems such as the ones treated in this work, the technique ought to be fine. Indeed, we have been successful in applying the method to spherocylindrical geometry at typical DNA parameters. These results will be presented in a future publication. One of the important aspects of this study has been to treat all the three different geometries on an equal footing and to calculate accurate electrostatic potential energy in each case using the same general principles. This study reveals that up to a certain limit, nonspherical macroions such as spherocylinder and spheroid shapes are also naturally overcharged like spherical macroions. In metastable ionized state, one can expect a long-range Coulomb attraction between two nonspherical macroions similar to the spherical macroion case. (52) Chodanowski, P.; Stoll, S. Macromolecules 2001, 34, 2320.

If the volume charge density is kept fixed, then for any macroion charge, the total electrostatic energy depends only on the geometry of the macroion. To observe the effect of only geometry change on electrostatic energy, the volumes of nonspherical macroions have been chosen to keep the same as the volume of a spherical macroion of radius R. This is also useful for comparison of electrostatic potential energies between macroions with different charges. We also have shown that the energy of neutral state is actually a continuous function of the shape factor τ. It is found that the total electrostatic energy of a spherical macroion with its neutralizing counterions (neutral state) can decrease to a certain minima (similar to gradual overcharging, first shown by Messina et al.13) by gradually deforming the macroion from its spherical shape, and in the particular case of spheroid (τ < 0), this energy decreases rapidly without turning. This clearly shows that not only by overcharging a spherical macroion but also by deforming it at a fixed charge can the total energy of the complex decrease. Finally, similar to fullerenes, the cell geometry of counterions is seen to be comprised of pentagons and hexagons, the only difference being that all pentagons and hexagons carry a counterion at their centers. Since this study establishes a technique with wider (spherical and nonspherical cases) applicability and basic information about the naturally overcharged nonspherical macroions, we expect to use this technique in the future to investigate areas related to screening such as proposed13,52 attraction between like-charged bundles of rodlike polyelectrolytes and determination of stable macroion clusters. It might be feasible to investigate whether there exists any stable ionic state for a couple of macroion in the presence of a salt. However, unlike MD, this technique can only depict the final stable state of the system so that studies of the intermediate metastable ionic states would be precluded. Acknowledgment. Partial support of an NSF Grant 0137273 is acknowledged. L.B.B. also wishes to acknowledge an internal institutional Grant from FIPI, University of Puerto Rico. A.K.M. expresses gratitude to Professors L. Fonseca and R. Messina for valuable discussions and e-mail correspondence, respectively, on some aspects of this work. Appendix A The “shape factor” τ is an important parameter in this study as it governs the geometry (shape) and size of the macroions. With the change of size (length and width), the volume of a spherocylinder or a spheroid is kept the same as that of a sphere of fixed radius R to maintain the volume charge density constant. For example, in the case of a spherocylinder, when the length L is changed, the radius r of the cylinder and also of the hemispherical endcaps is changed as shown below in eq 2A. We define the shape factor as τ ) L/r. Although τ can have any value, we have considered -1 e τ e 2 in this paper. In terms of τ, the aspect ratio of a sphercylinder is equal to L/(2r) ) τ/2. The radius of the cylinder and the end caps of the spherocylinder are equal. With the change of the length of the cylinder, this radius has to change to keep the volume the same as that of a sphere of radius R. This has been done in the following way:

πr2L + (4π/3)r3 ) (4π/3)R3

(1.A)

where r is the new adjusted radius and L is the length of

9612


Mukherjee et al.

the cylinder. With L ) rτ, eq 1A reduces to

r3τ + (4π/3)r3 ) (4π/3)R3 w r )

R ((3/4π)τ + 1)1/3 (2.A)

divide the charge in half and place a point charge Zm/2 at the locations “x” and also “-x” along the x-axis. If a counterion is located on the x-axis of the sphere, which we denote by the point X, it will experience an interaction distance relationship χX of magnitude

τ>0 Similarly, in the case of spheroids, we subtract a volume of a cross section (perpendicular to the z-axis) of length L from the sphere. The rest of the volume is now comprised of two equal sections of the original sphere whose radius can be adjusted to give a total volume equal to that of a sphere of fixed radius R. The adjusted radius is

r)

R [(3/4π)τ + 1 - (1/16)τ3]1/3

(3.A)

τ R, whose centers are L distance apart. Note that, if τ ) 0, either of the relations eq 2A or eq 3A, reduces to r ) R as expected for a sphere. Appendix B We show in this appendix that the internal distribution of the macroion charge directly affects the distribution of counterions on its surface. Consider as a reference a spherical macroion of radius R and charge Zm. If the charge is represented as a point at the center of the sphere designated by the point “O”, also the origin of the reference coordinate system, then the interaction energy with the counterion has the following distance relationship

χS ) Zm/R

(1.B)

where the subscript “S” denotes the spherical geometry. Clearly χS is valid for any point on the sphere. Let us now

χX )

(

)

(Zm/2) (Zm/2) Zm 1 + ) R+x R-x R 1 - δ2

(2.B)

where δ ) x/R. Clearly χX > χS, which means that the interaction of the counterion on the x-axis for the separated charge distribution is greater than that if the charge were located at the center of the sphere. Let us now use the same charge separation for the counterion located on the y-axis at the point which we denote by Y. The distance of the counterion from either charge at x or -x is now R/cos(φ) where the angle φ is defined by the points O, Y, and X, which is also expressed as x/R ) δ ) tan(φ). The interaction distance relationship in this case is

χY )

(Zm/2) (R/cos(φ))

+

(Zm/2) (R/cos(-φ))

)

Zm cos(φ) R

(3.B)

Clearly χY < χS, which means that the interaction of the counterion on the y-axis is less than that for the spherical shape. We can also use eqs 2.B and 3.B to deduce the effect of shape on the interaction energy for a noncentral macroion charge. In the metamorphosis from a sphere to a spherocylinder, we stretch the sphere along the x-axis and contract the y-dimension. Hence χX is now diminished in magnitude while χY increases in magnitude. In this geometric metamorphosis the counterion will summarily migrate from the point X to the point Y if the locations x and -x are fixed. LA0301590

Charge

Recommend Documents