Morphology of ion clusters in aqueous electrolytes

doi:10.1103/PhysRevE.77.031501

Phys Rev E Stat Nonlin Soft Matter Phys. Author manuscript; available in PMC 2008 July 15.

Published in final edited form as:

Phys Rev E Stat Nonlin Soft Matter Phys. 2008 March; 77(3 Pt 1): 031501.

Published online 2008 March 4. doi: 10.1103/PhysRevE.77.031501

PMCID: PMC2467436

NIHMSID: NIHMS51041

Morphology of ion clusters in aqueous electrolytes

Sergio A. Hassan

Sergio A. Hassan, Center for Molecular Modeling, DCB/CIT National Institutes of Health, Bethesda, Maryland 20892, USA;

See other articles in PMC that cite the published article.

Abstract

Formation of ion clusters in aqueous NaCl solutions at 25 °C is investigated with dynamics simulations in the 0.1–3M concentration range. The medium is found to be highly structured even at moderate concentrations, and clusters of over 20 ions are observed above ~2M. The medium can be viewed as a multicomponent fluid, composed of reacting particles with well-defined average populations, shapes, sizes, and charge distributions. Clusters show substantial morphological variations that are here characterized by their hydrodynamic radii and internal charge distributions. Clusters can be described as prolate spheroids in the subnanometer length scale. The hydrodynamic radius and the radius of gyration follow simple power laws with the number of ions in the cluster. Dipole moment distributions show characteristic peaks in the ~10–60 debye range that reflect conformational preferences modulated by electrostatic and liquid-structure forces.

I. INTRODUCTION

Electrolytes are commonly used in technological and industrial application, and are ubiquitous in biological systems. They are found in body fluids and tissues, and play a role in the highly regulated electrochemical balance in cells. Electrolytes affect physical and chemical properties of proteins and nucleic acids, the motility of individual cells (chemotaxis), the migration of simple multicellular organisms, and the survival of bacteria. These processes occur at a broad range of concentrations, from highly diluted to nearly saturated conditions.

Thermodynamic and transport properties of electrolyte solutions are well characterized [1]. However, their microscopic nature at ambient and physiological conditions is still poorly understood. Raman scattering experiments have suggested the presence of clusters in unsaturated, saturated, and supersaturated solutions [2,3]. More direct evidence of cluster formation has been obtained in recent years with dynamic light scattering techniques [4]. Assuming colloidal clusters of spherical shape, hydrodynamic radii were obtained for aqueous NaCl solutions at ~20 °C. Two well-defined size scales were identified: below ~1 nm and above ~100 nm [4].

The most direct evidence of cluster formation, however, is provided by computer simulations. They allow a detailed study of properties that cannot be scrutinized experimentally. Among these are the atomistic nature of cluster speciation, the properties of the liquid structure, and cluster morphology. Molecular dynamics (MD) simulations have been used to study structural and dynamic properties of concentrated solutions [5], the formation of small clusters in unsaturated aqueous electrolytes [6,7], early steps of salt nucleation in oversaturated solutions [8], and cluster formation in supercritical conditions [9].

Because of limitations in computer power earlier studies have been restricted to small systems and short simulation times. Recently, longer (10 ns) MD simulations provided direct evidence of clusters of up to six ions in 1m NaCl aqueous solutions at 25 °C [7] and permitted reliable calculation of thermodynamic data such as pair dissociation constants. In this paper, results from a set of 40-ns MD simulations of NaCl in water at 25 °C and 1 bar are reported. Salt concentrations are in the 0.1–3M range, which is typical in biological systems and related experimental studies. These relatively long simulations allow the observation of larger clusters, with 10–30 ions above ~2M, which start forming after ~10–20 ns. The simulations provide data for statistically reliable calculations, which permits identification of simple behaviors not yet reported in these systems.

II. METHODS

The simulations were performed with classical potentials in the canonical ensemble. A cubic cell of volume ~(46 Å)³ containing ~3500 nonpolarizable water molecules preequilibrated at 25 °C was used. The density was ≈0.993 g/cm³, consistent with a pressure of 0.1 MPa. Concentrations of 0.1M, 0.5M, 1M, 2M, and 3M were generated by randomly replacing water molecules with Na⁺ and Cl⁻ ions. Periodic boundary conditions were used (see additional details in [10]). The system was equilibrated for ~2 ns to allow ions to diffuse across the unit cell [11].

III. RESULTS

A cluster is defined as any array of ions such that (i) every ion is connected to at least one ion of opposite charge; two ions are said to be connected if they are separated by a distance smaller than 3.5 Å [12]; and (ii) every ion can be reached from any other ion through a path of consecutive connections. The following criterion applies to time evolution: If a cluster breaks into two fragments at time t but reforms before t+t, then the cluster experienced just a structural fluctuation, and did not break at t. In this study t is 2 ps [13]. Similarly, if two clusters coalesce at time t but separate again before t+t, then the clusters just collided, and did not react to form a single, larger cluster at t.

Statistical properties are derived from the cluster density ρ(n,q,t,c), defined as the number of clusters of size n (number of ions) and charge q per unit volume, at time t and concentration c. The clusters degree of formation is given by α(n,c) = c_n(c)/c, where equation M1

is the average molar concentration of clusters of size n; τ is the total simulation time, and N_A is Avogadro’s number. Figure 1 shows α_n vs n at different concentrations. The electrolyte is ~95% fully dissociated at 0.1M, and only ion pairs are present in the structured (associated) phase. The degree of dissociation decreases sharply as salt is added, as shown in the inset of Fig. 1 (n=1). At 3M the electrolyte is only ~45% dissociated, and the presence of very large clusters (n>20) can be observed as rare events. Although the largest clusters are observed infrequently and may have little effect on the system’s thermodynamics, they may play a role in nucleation and aggregation in longer time regimes or higher concentrations. The degree of formation of ion pairs increases as salt is added, but appears to reach a maximum around ~2M (Fig. 1 inset; n=2). Similar behavior is apparent from the data reported in [9], but here the maximum occurs at lower concentrations. The simulations show that large clusters can coexist in the relatively small space of the simulation cell. Above 2M, clusters of 10–20 ions are present simultaneously for several picoseconds before they break apart or grow further.

FIG. 1

Degree of formation α_n as a function of cluster size n (number of ions in the cluster) at different concentrations c (in moles/liter). The inset shows α_n as a function of c for dissociated ions (n=1) and for ion pairs (n=2).

The electrolyte undergoes significant structural changes in the course of the simulation. It fluctuates between highly structured and highly dissociated, with characteristic periods of ~0.5–1 ns. These changes can be quantified by the probability P_s that a fraction s of the electrolyte is in the structured phase. This probability can be defined as the time ratio P_s(s) = τ⁻¹∫G[s−s(t)]dt, where G(x) is 1 if x=0, and 0 otherwise; s(t) =1− f(t), and f(t) =c₁(t)/2c is the fraction of fully dissociated electrolyte at time t and concentration c. The calculation shows that the largest fraction s_M of structured electrolyte is ~35% at 0.1M, while the most probable fraction s₀ corresponds in this case to the lowest degree of structuredness, s_m, which is full dissociation. Fractions s_M, s₀, and s_m increase with salt, such that at 3M the electrolyte is highly structured at times, fluctuating between s_m ~45% and s_M ~70%, with a most probable fraction of s₀ ~55%. This shows that a strong electrolyte at room temperature can be highly associated even at moderate concentrations. It can be viewed as a multicomponent fluid, composed of reacting species with well-defined average populations. These particles vary in size, shape, and charge, as described below.

Cluster morphology shows substantial variation within classes (defined by n and q). Clusters are far from spherical, and present no crystal-like substructures. A measure of compactness can be defined as Ψ = γ/γ_c, where γ is the water accessible surface area (ASA) of a cluster in solution, and γ_c is the ASA of the most compact cluster that can be formed with the same ions. These compact clusters are taken as the smallest clusters seen in the NaCl crystal structure. This is the most efficient packing that would be physically expected in solution for this salt. Figure 2 shows γ vs n, which is found to be independent of concentration. γ is calculated numerically using the Lee-Richards definition [14] with a probe radius of 1.4 Å. It closely follows a power law of the form γ∞n^1/1.4 within the entire size range. Inset (a) of Fig. 2 shows Ψ vs n (solid, thick line), and suggests important deviations from ideal compactness (Ψ=1). To compare, Ψ is also given (solid, thin) for clusters with the lowest degree of compactness, taken here as a linear array of the ions in the crystal. Therefore, NaCl clusters in solution have ASA halfway between that of a cylinderlike and a spherelike array, and can then be viewed more simply as prolate spheroids. Yet, it is common practice to simplify the analysis of experimental data by invoking spherical colloids. The inset (b) of Fig. 2 shows effective radii calculated as R_γ=(γ/4π)^1/2 − 1.4 Å and R_c = (γ_c/4π)^1/2 − 1.4 Å. Thus, clusters in solutions have effective radii in the ~3–8 Å range and follow a power law R_γ∞n^1/d_γ, with dimension d_γ≈ 2.15. The simplest way to distort a sphere of radius R_c with ASA γ_c into a prolate spheroid of equatorial radius h_m and polar radius h_M with ASA γ, is to define h_m = R_c and h_M = (1+λ)R_c, where λ is a single parameter that quantifies the continuous deformation. The surface area of the spheroid is equation M2

, where the ellipcity ε is defined by ε ²=1−(h_m/h_M)². The problem is reduced to finding solutions of the equation 2(R_γ/R_c)²=1+(1+λ)sin⁻¹(ε)/ε for each value of n. The inset (b) of Fig. 2 shows λ ⁻¹ vs n (square symbols), and an exponential fitting λ ≈ 1.2–1.4 exp(−n/14) (dotted line). The ratio h_m/h_M decreases with n, so clusters sphericity decreases with size.

FIG. 2

Water accessible surface area γ as a function of cluster size n. Bars are standard deviations; statistical errors are small and omitted. The inset (a) shows the compactness parameter ψ (solid thick: simulation data; thin: linear arrange; (more ...)

Hydrodynamic radii R_H of colloidal particles are accessible experimentally using dynamic light scattering. Calculation of R_H by directly measuring mean square displacements is statistically unreliable due to the relatively short clusters lifetimes. Therefore, R_H is estimated here using the Bloomfield equation [15] for bodies of arbitrary shape, which is a generalization of the Kirkwood formula [16] for subunits of equal size. R_H of a cluster k of size n is then given by

where ρ_i is the hydrodynamic radius of ion (subunit) i, and r_ij is the distance between ions i and j; left angle bracket

_k denotes average over all conformations of cluster k. R_H is calculated as equation M4

, where N_n is the total number of unique clusters of size n observed in the course of the simulations. R_H vs n is plotted in Fig. 3, and found to be independent of concentration [17]. R_H follows closely a power law, R_H ∞n¹^/^d_H, within the entire cluster size range, with d_H ≈2.45. Plotted in the inset (a) of Fig. 3 is R_H vs R_γ; it shows little variations with size and a linear dependence with a ratio R_H/R_γ≈0.66.

FIG. 3

Hydrodynamic radius R_H as a function of cluster size. The inset (a) shows the linear dependence between R_H0 and R_γ. The inset (b) shows the linear dependence between R_H and the radius of gyration R_G. Bars are standard deviations; statistical errors (more ...)

Also accessible experimentally is the radius of gyration R_G, which can be measured using static light scattering. Inset (b) of Fig. 3 shows R_H vs R_G. Again the dependence is linear, with a ratio R_H/R_G ≈ 0.72 for n>5. The simulations do not show the formation of nucleating centers triggering irreversible aggregation. However, both the ratio R_H/R_G and the power-law behavior of R_H and R_G with n are similar to those observed in fractal aggregation of colloidal silica [18] and gold [19]. Independent measurement of R_H and R_G yielded R_H/R_G ≈ 0.72 in the limit of large silica aggregates, while Hausdorff dimensions of ~2.1 (silica) and ~2.0 (gold) have been estimated. The proportionality between R_H and R_G in fractal aggregates has been discussed in [20].

Cluster charge varies between −4e and +4e, depending on salt concentration, and will be studied elsewhere. The interest here is in higher moments of the cluster charge distribution, which displays substantial variations as well and provides information on the cluster morphology. Dipole moments are calculated here with respect to the clusters centers of mass, which can be far removed from the corresponding centers of charge (from which point μ=0). Thus, a cluster can be viewed as a liquid-excluding spheroidal body with charge q and dipole μ located at its center of mass. This interpretation may be useful to quantify clusters electrostatic effects on a distant solute (e.g., a protein or DNA molecule) through a multipole expansion of the potential. Bulk and nonbulk effects of concentrated electrolytes on biomolecular interactions have been studied in [21].

A distribution function p(μ|n,q) can be defined such that dP_n_,_q(μ)=p(μ|n,q)dμ is the conditional probability that a cluster of size n and charge q has a dipole moment of magnitude between μ and μ +dμ. Figure 4 illustrates the main features of p that highlight the preferred conformations in the case of small clusters. These conformations are determined by both electrostatic and liquid-structure forces. For n=2 the distribution reflects the variations of the interionic distance, yielding a peak at μ₂ ~12 D. For n=3, the additional internal degrees of freedom lead to broader distributions, with characteristic peaks in the μ₃ ~7–14 D range, depending on charge. This dependence stems from the difference in Na⁺ and Cl⁻ radii that affect the angle distribution at the central ion. At n=4, both extended and cyclic conformations are observed for neutral clusters. Extended clusters generate two peaks, at μ₄ ~12 D and μ₄ ~24 D; cyclic clusters yield one peak at a small value of μ₄ ~2.5 D. For n=4, charged clusters display less conformational variations because a central ion is surrounded symmetrically by three ions. This leads to relatively narrow distributions, and small dipoles in the μ₄ ~5–10 D range. As cluster size increases, the internal degrees of freedom lead to an increasing number of overlapping peaks. For n=5 two well-defined peaks appear in the μ₅ ~15–30 D range, depending on charge. Additional peaks can be inferred but are difficult to resolve visually. For n =6, neutral clusters show two well-resolved peaks in the μ₆ ~10–25 D range, but the distribution also suggests other preferred conformations. A decomposition of p as a sum of n Gaussians yields good fitting (not shown). This decomposition suggests that extended clusters shift the distribution towards larger values of μ, with a maximum peak at μ₆ ~40 D; partially cyclic clusters lead to smaller dipole moments, with a possible peak at μ₆ ~7 D. For n=10, a number of strongly overlapping peaks in the μ₁₀~20–80 D range are apparent that cannot be resolved. A general observation is that the locations of the peaks are independent of the salt concentration. The average dipole moment for each class is calculated as left angle bracket

μ_n_,_q

=∫ μ_n_,_qp(μ|n,q)dμ, and plotted in Fig. 5 for q=0 and q= ±1. It follows a linear behavior left angle bracket

μ_n_,_q

≈ a_q + b_qn, with a_q ≈ 5–10 D and b_q ≈ 2–3 D, and span the ~10–60 D range. Cluster survival times are in the ~10–100 ps range, depending on size and charge. Therefore, these large dipoles may have an important effect on the structure and dynamics of water, and profoundly affect the thermodynamics of the solution and of solvated biomolecules.

FIG. 4

Conditional probability distributions p as a function of clusters dipole moments μ (in debye) for small clusters. Distributions are normalized as ∫p(μ|n,q)dμ=1, and are shown to scale for each n. Charges are in units of (more ...)

FIG. 5

Average dipole moments left angle bracket

as a function of the cluster size. Bars are standard deviations; statistical errors are small and omitted. For n>25 (q=0) and n>20 (q= ±1) the trend deteriorates due to insufficient (more ...)

Acknowledgments

This study utilized the high-performance computer capabilities of the Biowulf PC/Linux cluster at the NIH. This work was supported by the NIH Intramural Research Program through the Center for Information Technology. The author thanks Peter Steinbach for reading the manuscript.

References

1. Robinson RA, Stokes RH. Electrolyte Solutions: The Measurement and Interpretation of Conductance, Chemical Potential and Diffusion in Solutions of Simple Electrolytes. Academic; New York: 1959.

2. Cerreta MK, Berglund KA. J Cryst Growth. 1987;84:577.

3. Rusli IT, Schrader GL, Larson MA. J Cryst Growth. 1989;97:345.

4. Georgalis Y, Kierzek AM, Saenger W. J Phys Chem B. 2000;104:3405.

5. Du H, Rasaiah JC, Miller JD. J Phys Chem B. 2007;111:209. [PubMed]

6. Degreve L, da Silva FL. J Chem Phys. 1999;111:5150.

7. Chen AA, Pappu RV. J Phys Chem B. 2007;111:6469. [PubMed]

8. Zahn D. Phys Rev Lett. 2004;92:040801. [PubMed]

9. Sherman DM, Collongs MD. Geochem Trans. 2002;3:102.

10. Simulations were carried out with the CHARMM program [22] using particle mesh Ewald summations with optimized parameters as reported in [23]. Electrostatic and van der Waals parameters of ions and the TIP3P water models are as implemented in the CHARMM force field. The water internal degrees of freedom were fixed with the SHAKE algorithm, and a time step of 2 fs was used. Data were collected every 1 ps for analysis.

11. Concentration-dependent self-diffusion coefficients D were calculated from the Einstein-Smoluchowski equation D= left angle bracket

|r(t)−r(t)|² right angle bracket

/6t at long t, where left angle bracket

denotes the average over all ions and over all time origins t₀ (0< t₀ <τ). The equilibration time was estimated from the diffusion of the slowest ions (Na⁺ at 3M, in which case D≈1.5 × 10⁻⁵ cm²/sec).

12. A rather restrictive definition of ion cluster is adopted here where l₀ ≈ 3.5 Å is the position (independent of concentration) of the minimum in the Na⁺-Cl⁻ pair correlation functions g(r) calculated from these simulations. A softer definition would use l₀ ≈ 6.0 Å, the position of the second minimum, thus allowing solvent-separated pairs within clusters.

13. The calculation of self-diffusion coefficients [11] implies that the slowest ions in bulk require ~2 ps to travel a distance equal to the radius of a water molecule (~1.4 Å). Thus, ~2 ps is taken here as the maximum time that an ion is allowed to move away from a cluster without exchanging position with a water molecule.

14. Lee B, Richards FM. J Mol Biol. 1971;55:379. [PubMed]

15. Bloomfield V, Dalton WO, van Holde KE. Biopolymers. 1967;5:135. [PubMed]

16. Kirkwood JG. J Polym Sci. 1954;12:1.

17. Hydrodynamic radii ρare calculated from the Stokes-Einstein relation ρ=kT/6πηD, where η is the viscosity coefficient of pure water. D’s were calculated as in [11]; η=0.8904 cp at 25 °C [24]. The calculations yielded ρ(Na⁺) ≈ 1.838 Å and ρ(Cl⁻) ≈ 1.206 Å.

18. Wiltzius P. Phys Rev Lett. 1987;58:710. [PubMed]

19. Weitz DA, et al. Phys Rev Lett. 1985;54:1416. [PubMed]

20. van Saarlos W. Physica A. 1987;147:280.

21. Hassan SA. J Phys Chem B. 2005;109:21989. [PMC free article] [PubMed]

22. Brooks BR, et al. J Comput Chem. 1983;4:187.

23. Hassan SA. J Phys Chem B. 2004;108:19501.

24. Weast RC. CRC Handbook of Chemistry and Physics. CRC Press; West Palm Beach: 1978.