Hybrid Monte Carlo and continuum modeling of electrolytes with concentration-induced dielectric variations

The distribution of ions near a charged surface is an important quantity in many biological and material processes, and has been therefore investigated intensively. However, few theoretical and simulation approaches have included the influence of concentration-induced variations in the local dielectric permittivity of an underlying electrolyte solution. Such local variations have long been observed and known to affect the properties of ionic solution in the bulk and around the charged surface. We propose a hybrid computational model that combines Monte Carlo simulations with continuum electrostatic modeling to investigate such properties. A key component in our hybrid model is a semianalytical formula for the ion-ion interaction energy in a dielectrically inhomogeneous environment. This formula is obtained by solving for the Green's function Poisson's equation with ionic-concentration-dependent dielectric permittivity using a harmonic interpolation method and spherical harmonic series. We also construct a self-consistent continuum model of electrostatics to describe the effect of ionic-concentration-dependent dielectric permittivity and the resulting self-energy contribution. With extensive numerical simulations, we verify the convergence of our hybrid simulation scheme, show the qualitatively different structures of ionic distribution due to the concentration-induced dielectric variations, and compare our simulation results with the self-consistent continuum model. In particular, we study the differences between weakly and strongly charged surfaces and multivalencies of counterions. Our hybrid simulations conform particularly the depletion of ionic concentrations near a charged surface and also capture the charge inversion. We discuss several issues and possible further improvement of our approach for simulations of large charged systems.

Figure 1

Schematic representation of a charged colloid immersed in an electrolyte. Cations and anions are represented by solid red and blue balls.

Figure 2

The autocorrelation function (ACF) with the surface charge density (a) $\sigma =-0.318e/{\mathrm{nm}}^2$ and (b) $\sigma =-1.214e/{\mathrm{nm}}^2$, and with different ionic valences.

Figure 3

(a) The integrated charge distribution functions (ICDFs) for three different rates of updating dielectric coefficient function: every 50, 1000, and $10000$ MC cycles. (b) The evolution of dielectric coefficient function at different iteration loops with the dielectric coefficient updated every 1000 MC cycles. Both (a) and (b) are calculated with the surface charge density $\sigma =-1.214e/{\mathrm{nm}}^2$ and with 3:1 electrolyte (trivalent counterions and monovalent coions).

Figure 4

Counterion and coion radial distribution functions (RDFs) of the distance to the surface of macroion. RDFs for monovalent, divalent, trivalent counterions, and coions with the ionic-concentration-dependent dielectric function $\varepsilon (c\left(r\right))$ are denoted by $g_{+},g_{2+},g_{3+},$ and $g_{-}$, respectively. RDFs for monovalent, divalent, trivalent counterions, and coions with a constant dielectric coefficient $\varepsilon =80$ are denoted by $g_{+}^0,g_{2+}^0,g_{3+}^0,$ and $g_{-}^0$, respectively. The anion concentration is fixed to be 100 mM. (a) 1:1 electrolyte and $\sigma =-0.318e/{\mathrm{nm}}^2$. (b) 1:1 electrolyte and $\sigma =-1.214e/{\mathrm{nm}}^2$. (c) 2:1 electrolyte and $\sigma =-0.318e/{\mathrm{nm}}^2$. (d) 2:1 electrolyte and $\sigma =-1.214e/{\mathrm{nm}}^2$. (e) 3:1 electrolyte and $\sigma =-0.318e/{\mathrm{nm}}^2$. (f) 3:1 electrolyte and $\sigma =-1.214e/{\mathrm{nm}}^2$.

Figure 5

The comparison of caion (anion) density profiles by MC simulations and SCCM model for a 1:1 dilute electrolyte (13 mM) around a neutral $(\sigma =0)$ macroparticle, where $g_i^0$ is $g_{+}^0$ or $g_{-}^0.$

Figure 6

Counterion and coion radial distribution functions (RDFs) of the distance to the macroion surface from three different models: MC simulations ($g_i$), self-consistent continuum model ($g_i^{\text{CM}}$), and PB model ($g_i^{\text{PB}}$). Two cases are considered here: (a) a $2:1$ electrolyte with a low surface charge density $\sigma =-0.318e/{\text{nm}}^2$ (a); (b) a $1:1$ electrolyte with a high surface charge density $\sigma =-1.214e/{\text{nm}}^2$.