The solvation and ion condensation properties for sulfonated polyelectrolytes in different solvents—a computational study

The evaluation of statistical mechanics methods on the solvent and ion distribution function allows important insights into the preferential binding behavior in terms of the corresponding Kirkwood–Buff theory, which was introduced in the early 1950s [35, 36]. The radial distribution function of molecules or atoms β around solutes α can be expressed by

$$\begin{equation} g_{\alpha\beta}(r) = \frac{\rho_{\beta}(r)}{\rho_{\beta,\infty}}, \end{equation} \tag{ 1 }$$

where ρ_β(r) denotes the local density of β at a distance r around the solute, and ρ_β,∞ is the global density in the bulk phase [37]. The integrated radial distribution function gives the cumulative number distribution function

$$\begin{equation} f_{\alpha\beta}(r) = 4\pi\int_0^r g_{\alpha\beta}(r)\,\mathrm{d}r, \end{equation} \tag{ 2 }$$

which yields an estimate for the number of molecules of type β within a radius r around molecule α [38]. The Kirkwood–Buff integral is given by the integration of equation (1)

$$\begin{equation} G_{\alpha\beta} = \lim_{R\rightarrow\infty} G_{\alpha\beta}(R) = \lim_{R\rightarrow\infty}\int_{r=0}^{r=R}4\pi r^2(g_{\alpha\beta}(r)-1)\;\mathrm{d}r, \end{equation} \tag{ 3 }$$

where the above relation is valid in the limit of R = ∞ [35, 36, 39–41]. Due to the presence of a finite box length in computer simulations, one typically defines the radial Kirkwood–Buff parameter as

$$\begin{equation} G_{\alpha\beta}(r) = \int_{0}^{r}4\pi r^2(g_{\alpha\beta}(r)-1)\;\mathrm{d}r, \end{equation} \tag{ 4 }$$

where the lower integration limit defines the molecular surface and the upper value is given at the point where the values of G_αβ(r) converge [39, 40]. Equation (3) can be used to calculate the excess coordination number of atoms β (sodium ions) around α (polyelectrolyte) via [35, 36, 39–41]

$$\begin{equation} N_{\beta}^{xs} = \rho_{\beta,\infty}G_{\alpha\beta}= \rho_{\beta,\infty}\lim_{R\rightarrow\infty} G_{\alpha\beta}(R), \end{equation} \tag{ 5 }$$

which allows us to evaluate the preferential binding coefficient ν_βγ(R). Under the assumption of finite distances r, the radial preferential binding coefficient is given by

$$\begin{equation} \nu_{\beta\gamma}(r) = \rho_{\beta,\infty}(G_{\alpha\beta}(r)-G_{\alpha\gamma}(r)) = N_{\beta}^{xs}(r) - \frac{\rho_{\beta,\infty}}{\rho_{\gamma,\infty}}N_{\gamma}^{xs}(r), \end{equation} \tag{ 6 }$$

where the index γ represents solvent molecules (water, chloroform or DMSO), while α and β denote the polyelectrolyte and the sodium ions. A positive value for the preferential binding coefficient of equation (6) implies an energetically favorable binding of sodium ions [39, 40] due to

$$\begin{equation} \Delta F_{\beta\gamma}^{\mathrm{t}}(r) = -RT \nu_{\beta\gamma}(r), \end{equation} \tag{ 7 }$$

which illustrates the connection of the preferential binding coefficient with the transfer free energy ΔF^t_βγ multiplied with the molar gas constant R and temperature T. The transfer free energy gives an estimate for the amount of energy that is needed to transfer a sodium ion from a point r to the close vicinity of the polyelectrolyte surface.

Several methods have been proposed over the years ([37, 43, 44] and references therein) to obtain an estimate for the binding free energy between two species. It has been often discussed that the estimation of free energy differences is a challenging task in computer simulations. A well-defined sampling of the configurational space is mandatory. An easier option that is applicable for solvent–solute interactions is the estimation of the binding free energy F(r) at the distance r between two species via the pair radial distribution function g_αβ [37, 45]

$$\begin{equation} \Delta{F}(r) = -k_{\rm B}T\,\log\,(g_{\alpha\beta}(r))+C, \end{equation} \tag{ 8 }$$

where the constant C is given by reference energy values. With the binding free energies, the solvation properties for different solvents can be estimated and solvophobic as well as solvophilic properties can be identified.

To study the influence of different solvents on the conformational properties of the polyelectrolyte, we have determined the radius of gyration, the end-to-end distance and the solvent accessible surface area for mod-s220 Na⁺. The radius of gyration R_g is given by

$$\begin{equation} R_{\mathrm{g}}^2 = \frac{1}{2n^2}\sum_{i,j}^n(\vec{R}_i-\vec{R}_j)^2, \end{equation} \tag{ 9 }$$

where n denotes the number of monomers and $\vec {R}_i,\vec {R}_j$ the positions for different monomers i,j [56]. Another expression that can be used to estimate the conformational behavior is given by the end-to-end distance

$$\begin{equation} R_{\mathrm{e}}^2 = (\vec{R}_1-\vec{R}_n)^2, \end{equation} \tag{ 10 }$$

which connects the position of the first monomer $\vec {R}_1$ to the last monomer $\vec {R}_n$ [56]. The values for the different solvents are shown in table 1. It can be clearly seen that the average values for R_g and R_e are comparable for DMSO and water. Nevertheless, the values for water are more fluctuating in regards to a more-pronounced standard deviation. One reason for this behavior is the stronger fluctuation of the number of condensed counterions in the presence of water compared to DMSO and chloroform, as will be discussed in the next subsection. In comparison to the values for DMSO and water, it can be concluded that chloroform is a significantly poorer solvent for the polyelectrolyte. The reason for the significant increase of the polyelectrolyte size in DMSO and water can be related to favorable solvent–polyelectrolyte interactions in agreement with the Flory–Huggins theory [56].

The polarity of the solvent in terms of different dielectric constants and their influence on the polyelectrolyte conformation can be also investigated by the amount of the exposed polar solvent accessible surface area Σ₊, which is calculated by taking into account the solvent accessible surface areas of the sulfur and oxygen atoms. It can be seen in table 1 that Σ₊ becomes larger in the order Σ₊(water) >Σ₊(DMSO) >Σ₊(chloroform). The coincidence between the polarity of the solvent as given by the dielectric constant and the amount of the polar solvent accessible surface area is obvious. In contrast, no trend can be recognized for the behavior of the apolar solvent accessible surface area Σ₋ (all atoms except sulfur and oxygen) with respect to the polarity of the solvent. The order is given by Σ₋(DMSO) > Σ₋(chloroform) >Σ₋(water). This behavior indicates that the larger solvent accessible surface area of the polar groups in response to the polarity of the solvent is the main driving factor for extended configurations. Our results further evidence that the interactions of the polar groups with water are the most favorable ones, as can be seen by the largest polar solvent accessible surface area. In contrast, the large polarity of water avoids an interaction with the apolar groups, which explains the smaller values for the apolar solvent accessible surface area due to the hydrophobic effect [20]. Hence, the polyelectrolyte has to find a balance between the favorable and unfavorable exposed solvophilic and solvophobic solvent accessible surface area. We assume that the intermediate polarity of DMSO between water and chloroform, as given by the dielectric constant  ≈ 47, favors a balanced increase in the polyelectrolyte size for the polar as well as the apolar groups, which results in the largest total solvent accessible surface area Σ_t compared to the other solvents.

The ratio of the apolar to the total solvent accessible surface area R = Σ₋/Σ_t is shown in figure 3. It can be easily seen that the ratio R is significantly dependent on the polarity of the solvent. We have found values of R(water) = 0.40 ± 0.01, R(DMSO) = 0.44 ± 0.01 and R(chloroform) = 0.52 ± 0.02. Thus, a decreased ratio of the apolar regions can be observed for polar solvents, which validates the presence of the solvophobic or hydrophobic effect. The opposite trend is also valid for the ratio of the polar to the total solvent accessible surface area by R₋ = 1 − R. It becomes evident that the increase in the polyelectrolyte size is given by a subtle interplay between the apolar and polar regions of the polyelectrolyte and the corresponding response to the polarity of the solution. These findings are in good agreement with protein folding theory and the aggregation of hydrophobic particles, where it was assumed that the minimization of hydrophobic regions in terms of the hydrophobic collapse is the initial starting point for the formation of protein secondary and tertiary structures as well as self-aggregation of hydrophobic and amphiphilic molecules [19, 20].

Zoom In Zoom Out Reset image size

Another important factor that influences the polyelectrolyte conformation is given by the number of condensed counterions. A strong condensation behavior corresponds to an efficient screening of electrostatic interactions between the charged groups, which results in weak repulsive forces between the monomeric charges [15, 61]. In regards to the qualitative results in figure 2, it can be seen that the largest number of condensed counterions occurs for chloroform followed by water. Obviously, the presence of DMSO leads to the smallest number of condensed sodium ions. Thus, the slightly more pronounced larger polyelectrolyte conformation size for DMSO, as given by the values of table 1, can be also related to stronger repulsive electrostatic interactions between the likewise charged $-\!\!\!-{\rm SO}_{3}-\!\!\!-$ groups due to the smaller number of screened charges and the lower dielectric constant of the solvent. We will evaluate this point in the next section in more detail by the introduction of a quantitative analysis method.

The binding free energy between solvent molecules and the different polar and apolar groups is shown in figure 4. It can be clearly seen that DMSO leads to the strongest binding behavior in the first solvation shell around the $-\!\!\!-{\rm SO}_{3}-\!\!\!-$ (SO3) group. The well-pronounced energetic minimum is located at 0.4 nm at the first solvation shell with an energy barrier of 1.5 kJ mol⁻¹, which prevents the release of solvent molecules to the second solvation shell at 0.8 nm. The binding energy compared to the bulk solution at 1.5 nm is negative with ΔF = −0.4 kJ mol⁻¹, which indicates favorable interactions. The behavior for water is slightly different. The first energy minimum occurs at 0.3 nm with ΔF ≈ −0.1 kJ mol⁻¹ and an energy barrier of roughly 0.8 kJ mol⁻¹. At this distance, most of the hydrogen bonds between the polyelectrolyte and water are present. The various small peaks within a distance of 0.3–0.45 nm can be related to different conformations of the water molecules around the SO3 group. In summary, it can be concluded that the overall binding energy of DMSO to the SO3 group is larger that of water. As was expected, the binding energy of chloroform is unfavorable for all distances due to the apolar properties of the solvent and the polar characteristics of the SO3 group. A nearly identical behavior to the SO3 group can be also observed for the SO2 group with less pronounced energy barriers and energy minima.

Zoom In Zoom Out Reset image size

In addition, we have also analyzed the binding properties of the different solvents to the phenyl ring (carbon). It can be seen that the binding free energy is positive for all solvents. However, the most pronounced free energy minimum can be observed for chloroform at a distance of 0.6 nm. Hence, although the values are positive, the weakly pronounced affinity of association between apolar solvents and apolar groups becomes visible. In regards to the above-discussed results, it can be concluded that the molecular details of the polyelectrolyte solvent contact include a complex solvation mechanism due to alternating polar and apolar groups in the polyelectrolyte backbone. However, all estimated binding free energies are in good agreement with the above-discussed conformational behavior.

To further evaluate the binding properties of DMSO to the polyelectrolyte, we have calculated the spatial distribution of the solvent's center of mass for DMSO around an arbitrary middle monomer of the polyelectrolyte, as shown in figure 5. It can be clearly seen that a large fraction of DMSO molecules accumulates between the three sulfur-containing groups SO3 and SO2. It is further evident that DMSO molecules aggregate around the carbon ring and the SO3 group. Therefore, we assume that electrostatic interactions such as dipole–dipole interactions between DMSO and the partially charged groups of the polyelectrolyte are the main driving factor for the solvation properties of DMSO. For the investigation of the DMSO orientation, we have evaluated the radial distribution function g_S−C/O between all sulfur atoms of the SO3 and SO2 groups and the corresponding carbon and oxygen atoms, respectively, of DMSO. The results are presented on the right side of figure 5. It can be clearly seen that a well-ordered orientation is given for the methyl groups (carbon atoms) of DMSO toward the sulfur atoms of the polyelectrolyte with the highest peak at 0.45 nm. It can be concluded that the DMSO molecules and specifically the slightly positively charged methyl groups and the central sulfur atom are electrostatically attracted by the oxygen atoms of the SO3 and SO2 groups. It has to be noted that the corresponding orientation analysis for water is not that well pronounced due to the larger number of condensed sodium ions that electrostatically interact with the solvent molecules. We will come back to this point in more detail in the next section.

Zoom In Zoom Out Reset image size

We have further estimated the cumulative number distribution function of solvent molecules around the polyelectrolyte surface area according to equation (2), and have multiplied it with the molecular volume of the distinct solvent molecules, as shown in figure 6. The multiplication by the molecular volume allows us to take excluded-volume effects into account. We have determined the molecular volumes of V_m(water) ≈ 0.065 nm⁻³, V_m(DMSO) ≈ 0.217 nm⁻³ and V_m(chloroform) ≈ 0.232 nm⁻³. It can be clearly seen that DMSO and water occupy a larger volume around the polyelectrolyte compared to chloroform, which can be related to favorable interactions. After the first solvation shell for water and DMSO around 0.4 nm, it can be seen that DMSO and water have a nearly identical occupied volume. The main differences occur due to the ordering behavior in the second solvent shell of DMSO, as given for values around 0.6 nm. A strong depletion effect can be therefore observed for DMSO within a distance of 0.5–0.6 nm, which corresponds to the broad and pronounced transition and energy barrier, respectively, of the first to the second solvation shell in agreement with the results of figure 4.

Zoom In Zoom Out Reset image size

As we have shown qualitatively in figure 2, the ionic condensation behavior increases for the case where DMSO < water < chloroform. This observation is validated by figure 7, where the cumulative number distribution function of sodium ions around the polyelectrolyte surface according to equation (2) is shown. We have identified the fraction of condensed sodium ions by considering the number of ions within the corresponding Bjerrum length. A nearly full association of sodium ions around mod-s220 Na⁺ can be observed in the presence of chloroform. Due to the apolar character of chloroform, the counterions tend to associate with the charged polyelectrolyte groups. In addition, the very large Bjerrum length of λ_B ≈ 25 nm, due to the low dielectric constant in the chloroform solution, leads to strong association properties. The observed behavior for chloroform can be explained by the Manning–Oosawa counterion condensation theory [58–60, 62]. It was originally proposed that counterion condensation occurs if Γ = λ_B/l_c ≫ 1 where l_c is the distance between two neighboring charged polyelectrolyte groups. For chloroform, it is clearly evident that Γ ≫ 1, which results in the presented behavior.

Zoom In Zoom Out Reset image size

The situation is less clear if the condensation behavior for DMSO and water is considered. Figure 7 significantly illustrates that more ions are condensed in the presence of water compared to DMSO, although Γ_water < Γ_DMSO due to a smaller water Bjerrum length. If we calculate the fraction of condensed counterions by dividing the number of counterions within the distance of the Bjerrum length N_λB with the total number of counterions N in terms of r_c = N_λB/N, we derive a condensation constant of r_c ≈ 0.86 for water, while DMSO leads to a constant of r_c ≈ 0.38 (chloroform r_c = 1). The low condensation constant for DMSO may also explain the larger polyelectrolyte conformation size compared to water due to unscreened repulsive electrostatic interactions between the SO3 groups, as was discussed above. Furthermore, it has to be remarked that the number of condensed counterions in water is strongly fluctuating, as has been observed in the larger standard deviations compared to DMSO and chloroform. As a consequence, we can relate this effect to the larger fluctuations of the radius of gyration and end-to-end distance as were shown in table 1.

A main driving force for the ionic condensation behavior is given by the counterion binding free energy, which can be also expressed by the transfer free energy according to equation (7). The affinity for association at each distance can be calculated by the radial preferential binding coefficient ν_βγ(r) according to equation (6). The corresponding radial distribution functions have been evaluated by taking into account the sulfur atom of the SO3 group and the corresponding centers of masses for the solvent molecules and sodium ions, respectively. The values for the radial preferential binding coefficient are shown in figure 8. It can be clearly seen that the preferential binding coefficient is positive for all distances and solvents. This indicates that the transfer of ions from the bulk solution to the polyelectrolyte surface according to equation (7) is a favorable process. It becomes obvious that the transfer free energy for chloroform is roughly 4 times and 2.5 times, respectively, more favorable compared to DMSO and water. These large values clearly indicate that free counterions in chloroform are nearly absent in agreement with the above-discussed results. It can be additionally seen that the association of sodium ions around the polyelectrolyte in an aqueous solution is more favorable compared to DMSO, which is in agreement with the determined higher condensation constant r_c for water. This behavior has been also experimentally validated [63], where it was found that the transfer free energy of a sodium ion from an aqueous solution to DMSO is negative. With regards to these results, we conclude that the pronounced deviations in the ionic solvation properties for DMSO and water are responsible for the different observed ionic condensation constants.

Zoom In Zoom Out Reset image size

The detailed solvation behavior for the sodium ions can be evaluated in terms of the pair radial distribution function, as shown in figure 9. We have considered all sodium ions within a distance of 0.7 nm of the macromolecular surface, as condensed to the polyelectrolyte. It can be clearly seen that the number of water molecules in the first hydration shell, up to a distance of 0.32 nm, differs for condensed and free sodium ions. This can be mainly explained by the close distance to the polyelectrolyte in terms of steric interactions. We have further evaluated the cumulative number distribution function by the corresponding evaluation of the pair radial distribution functions according to equation (2), to determine the number of solvent molecules in the first shell. The corresponding results at 0.32 nm are given by N_c = 5.45 ± 0.02 water molecules for free sodium ions and N_c = 4.53 ± 0.14 water molecules for condensed sodium ions. The results for the free ions are in good agreement to a recent publication where this number has been determined to be around N_c = 5.7–5.8 [64]. Furthermore we have also evaluated the corresponding binding free energies of sodium ions to solvent molecules by using equation (8). Compared to the bulk solution at r ⩾ 1 nm where the radial distribution function converges, we have determined ΔF(water,free) ≈ −4.2 kJ mol⁻¹ for the free ions and ΔF(water,condensed) ≈ −3.3 kJ mol⁻¹ for the condensed ions. The corresponding analysis for DMSO (the middle sketch of figure 9) shows that there is really a negligible difference in the solvation behavior for condensed and free sodium ions in DMSO. It can be assumed that a full solvation shell is always around the sodium ions, which indicates the good solubility of sodium in DMSO. This is also clearly evidenced by the large binding free energies of condensed and free ions with ΔF(DMSO,free/condensed) ≈ −5.3 kJ mol⁻¹, which is around 1 kJ mol⁻¹ more favorable compared to water. If the number of DMSO molecules in the first solvation shell, up to a distance of 0.42 nm, is analyzed, we find nearly identical values for condensed and free ions with N_c = 5.97 ± 0.02 and 5.99 ± 0.02 DMSO molecules. These values are in qualitative agreement with literature results, which have been estimated by ab initio methods [65]. As a reason for the good solubility, it has been emphasized that the solvation properties of sodium ions in DMSO are extremely well pronounced. DMSO can be considered as a coordinating solvent in terms of significant electron pair donor abilities [66]. These properties can be also expressed by the so-called donor numbers [66, 67], which give an estimate of the solubility of an ion in a specific solvent. The donor numbers for water and DMSO are given by DN_H2O = 18 kcal mol⁻¹ and DN_DMSO = 29.8 kcal mol⁻¹, whereas chloroform shows a vanishing donor number [67] in agreement with our findings with regard to the condensation properties.

Zoom In Zoom Out Reset image size

It is evident that there is a strong difference in the solvation behavior for condensed and free sodium ions in chloroform, despite the low statistical accuracy due to the small number of free counterions at distances r ⩾ 0.6 nm. We have identified values of N_c = 5.63 ± 0.18 chloroform molecules for free, and N_c = 3.21 ± 0.14 chloroform molecules for condensed sodium ions. The binding energy is smallest compared to DMSO and water, and has been estimated to ΔF(chloroform,free) ≈ −1.1 kJ mol⁻¹ for the free ions and ΔF(water,condensed) ≈ −0.2 kJ mol⁻¹ for the condensed ions. Hence the small values for the binding free energy clearly indicate the poor solubility of sodium ions in a chloroform solution, which may also result in the large condensation constant for this solvent and which is in agreement with the donor number concept [67].

In this subsection, we present the dynamic behavior of the sodium ions in different solvents. Due to the fact that a nearly complete condensation of sodium ions in a chloroform solution has been validated, we purely restrict our analysis to water and DMSO. The sodium and solvent diffusion constants have been calculated via the mean-squared displacement [37, 42]

$$\begin{equation} D = \lim_{t\rightarrow\infty}\frac{\langle (\vec{R}(t)-\vec{R}(t_0))^2\rangle }{6t}, \end{equation} \tag{ 11 }$$

where $\vec {R}$ denotes the position of the sodium ions and the center-of-mass position for a solvent molecule at time t and t₀, respectively. The corresponding values are presented in table 2. It has to be noted that the calculated diffusion constant for the SPC/E water model is in excellent agreement with a previous publication [68]. The diffusion constants for DMSO are also in good agreement with additional results for united atom models [69, 70]. The comparison of the sodium diffusion constants between DMSO and water reveals that the value for water is roughly two times larger, in contrast to the pronounced difference between the solvent diffusion constants, which can also be explained by the larger condensation constant of sodium ions in an aqueous solution. In regards to a recent publication [5], it was proposed that the ionic dissociation constant can also be estimated by the ratio of the ion to the solvent diffusion constant via r_diss = D_Na⁺/D_s. For the validity of this relation, it has to be assumed that the diffusion constant for solvated and dissociated ions roughly agrees with the solvent diffusion constant. A dissociation constant of r_diss ≈ 0.15–0.2, which corresponds to a condensation constant of r_c,e = 0.8–0.85 for lithium ions in DMSO, has been estimated by the experiments for λ values around 8 [5]. As discussed in the last section, we have found a value for the sodium condensation constant of r_c,t ≈ 0.38, which is significantly smaller than the experimental prediction. Following the approach presented in [5], we derive a dissociation constant for sodium ions in DMSO of D_Na⁺/D_s = 0.31, which corresponds to a fictive condensation constant of 0.69. The difference between the differently derived condensation constants illustrates the fact that counterion condensation properties are hard to estimate by taking the diffusion constants into account.

An important estimator for the successful applicability of polyelectrolytes like mod-s220 Na⁺ in technological products is given by the sodium conductivity σ. Although our simulated system is not directly comparable to experimental setups where a much higher counterion concentration is envisaged, an estimate for the ion conduction quality can be achieved by analyzing the simulation results. A reasonable option for determining the sodium conductivity in experiments and also for computer simulations is the usage of the Nernst–Einstein relation for cations [71]

$$\begin{equation} \sigma_{\rm Na^+} = \frac{Nq^2}{Vk_{\rm B}T}D_{\rm Na^+}, \end{equation} \tag{ 12 }$$

where N is the number of sodium ions and q = 1e the corresponding charge. The corresponding values are shown in table 2. Surprisingly, it can be seen that the sodium conductivity in water is roughly six times higher than DMSO, although the counterion condensation constant in water is larger. Therefore it can be concluded that a good choice for a solvent to optimize ionic and ion conductivities is given by the realization of large diffusion constants with a small ionic condensation constant.