Discrete ion stochastic continuum overdamped solvent algorithm for modeling electrolytes

In this paper we develop a methodology for the mesoscale simulation of strong electrolytes. The methodology is an extension of the fluctuating immersed-boundary approach that treats a solute as discrete Lagrangian particles that interact with Eulerian hydrodynamic and electrostatic fields. In both algorithms the immersed-boundary method of Peskin is used for particle-field coupling. Hydrodynamic interactions are taken to be overdamped, with thermal noise incorporated using the fluctuating Stokes equation, including a “dry diffusion” Brownian motion to account for scales not resolved by the coarse-grained model of the solvent. Long-range electrostatic interactions are computed by solving the Poisson equation, with short-range corrections included using an immersed-boundary variant of the classical particle-particle particle-mesh technique. Also included is a short-range repulsive force based on the Weeks-Chandler-Andersen potential. This methodology is validated by comparison to Debye-Hückel theory for ion-ion pair correlation functions, and Debye-Hückel-Onsager theory for conductivity, including the Wien effect for strong electric fields. In each case, good agreement is observed, provided that hydrodynamic interactions at the typical ion-ion separation are resolved by the fluid grid.

Figure 1

Illustration of short-range potential and repulsive force due to Eq. (42). For ease of viewing, here we have set ${\tilde{x}}_m=0.7\sigma$. Note that the actual value used in the simulations, discussed in Sec. 5, is ${\tilde{x}}_m=0.25\sigma$.

Figure 2

Left: a 2D illustration of the MAC [54] discretization used for Eq. (50). Right: 2D illustration of the ${\tilde{\mathbf{\nabla }}}_h·\mathit{Z}$ term in Eq. (50a), which is the discretization of the term ${\mathbf{\nabla }}_r·\mathcal{Z}$ from Eq. (11). Random numbers $\mathit{Z}$ are generated at the faces of control volumes around the points where velocities $\mathit{v}$ are defined so that the divergence of the field $\mathcal{Z}$ can be calculated.

Figure 3

Radial distribution for a molarity of 0.1 M. The right plot shows the pair correlation for ions of like charge; also indicated is the Debye length. The left plot shows the pair correlation for ions of opposite charge. The solid line shows the approximate analytical solution given by Eq. (57), and the squares show the numerical result from DISCOS. These results were collected using a bin width of 0.025 nm, with a sample size of $50000$; the simulation was run for $100000$ steps, with the first $50000$ steps used for equilibration. This sample size resulted in negligible statistical error. For visual clarity, only every fourth point is displayed. The grid size was set such that the total diffusion of species $A$ was 14.7% wet and 85.3% dry, and for species $B$, 12.9% wet and 87.1% dry. This corresponds to $\mathrm{\Delta }r=1\text{nm}$, with a $32\times 32\times 32$ cell periodic domain, and 3946 ions. The time step was $\mathrm{\Delta }t=0.1$ ps.

Figure 4

Radial distribution for a molarity of 1.0 M. The right plot shows the pair correlation for ions of like charge; also indicated is the Debye length. The left plot shows the pair correlation for ions of opposite charge. The solid line shows the approximate analytical solution given by Eq. (57), and the boxes and circles show numerical result from DISCOS. The numerical parameters are the same as those described for Fig. 3, except that two different grid resolutions are shown. The black circles show the case where the grid size was set such that the total diffusion of species $A$ was 14.7% wet and 85.3% dry, and for species $B$, 12.9% wet and 87.1% dry. This corresponds to $\mathrm{\Delta }r=1\text{nm}$, with a $16\times 16\times 16$ cell periodic domain, and 4932 ions. The red squares show the case where the grid size was set such that the total diffusion of species $A$ was 58.9% wet and 41.1% dry, and for species $B$, 51.8% wet and 49.2% dry. This corresponds to $\mathrm{\Delta }r=0.25\text{nm}$, with a $32\times 32\times 32$ cell periodic domain, and 618 ions.

Figure 5

Wien corrections to the DHO relaxation effect and electrophoretic contribution as a function of $\mathfrak{a}=\left(q{\lambda }_{D}E\right)/\left(k_{\text{B}}T\right)$.

Figure 6

Relative difference in conductivity from DHO theory including Wien corrections for the strong-field case, as a function of wetness percentage, for a range of molarities; based on data from Table 1.

Figure 7

Left: mobility of particle 1 when a force is applied on particle 2, in the direction parallel to the vector connecting particles 1 and 2. Right: mobility of particle 1 when a force is applied on particle 2, in the direction perpendicular to the vector connecting particles 1 and 2. The red circles are values measured from a 100% wet simulation, the black crosses from a 50% wet simulation. In each case the values are normalized by the particle self-mobility, and the separation is shown in terms of the average particle radius of species $A$ and $B$, i.e., ${\overline{a}}_t=(a_t^A+a_t^B)/2, {\overline{a}}_w=(a_w^A+a_w^B)/2, {\overline{a}}_d=(a_d^A+a_d^B)/2$. Each set of data is compared to the result from the RPY tensor calculated using $a_w$, shown with the red and black solid lines. Also shown for comparison is the RPY result corresponding to the 75% wet case (dashed line). The approach of using the average particle radius has been used for simplicity because the two species are of similar size; we note that a polydisperse version of the RPY tensor is described in Ref. [67]. Simple fits of the pair mobility for the six-point Peskin kernel is given in Ref. [68].

Figure 8

Particle velocity due to electrostatic and hydrodynamic interactions between two ions, as a function of separation. The velocity $V_1$ is calculated using Eqs. (72, 73, 74), using a radius averaged from species $A$ and $B$, i.e., ${\overline{a}}_t=(a_t^A+a_t^B)/2, {\overline{a}}_w=(a_w^A+a_w^B)/2, {\overline{a}}_d=(a_d^A+a_d^B)/2$. A range of wet percentages are shown, the corresponding particle radii can be found with Eqs. (29) and (30) and the radii quoted at the start of Sec. 5. The measured average closest ion separations corresponding to concentrations 0.1 and 2 M are indicated with the black and red vertical lines.

Figure 9

Normalized particle velocity due to Eq. (74) (lines), and normalized conductivity simulation measurements from Table 1 (symbols), as a function of wet and dry ratio of species $A$. All values are normalized by the 0% wet value. Black signifies 0.1 M, red 2 M. The curves indicate the particle velocity, calculated at the average nearest particle separation measured from simulation (solid), and corresponding to random (dashed), and cubic lattice (dotted) configurations. Squares are zero-field conductivity measurements using the Einstein-Helfand formula [Eq. (70)], diamonds indicate weak-field conductivity measurements using the center-of-charge displacement [Eq. (68)], and crosses are the prediction of DHO theory [Eq. (61)].

Figure 10

Local correction function $F^{\mathrm{LC}}(x=r/\mathrm{\Delta }r)$ computed from the radial component of the near-field correction to the electrostatic force ${\mathit{F}}_{ij}^{\mathrm{LC}}$ for the Peskin four-point kernel (left panel) and the Peskin six-point kernel (right panel). Red circles are samples for different positions and orientations of the pair of ions, while dashed blue lines indicate the continuum theoretical near-field correction $F_{\text{Gaussian}}^{\mathrm{LC}}$ for Gaussian charges of standard deviation $\sigma$ indicated in the legend [see Eq. (A2)]. The black solid line in the left panel indicates the empirical mean used to tabulate $F^{\mathrm{LC}}$ in our method. The right panel includes results (black squares) for a spectral Poisson solver with a six-point exponential of a semicircle (ES) kernel [79] with parameter $\beta =1.8\times 6$.

Figure 11

Percent error due to loss of translational and rotational invariance in the immersed-boundary P3M method based on data shown in Fig. 10. In addition to the method used in this work, we include for comparison the improvement gained by using a fourth-order isotropic Poisson solver, as well as a spectral Poisson solver with a six-point ES kernel. (Left panel) Percent error in the radial component of the total electrostatic force due to loss of invariance. (Right panel) Magnitude of the nonradial component of the IB force $F^{\mathrm{P}}$ expressed as percent of the radial component. The symbols show the empirical mean over samples of pairs of points, while the dashed lines show the maximum value over the samples.

Figure 12

Left: comparison of cation molarity between DISCOS (red) and MD (black), for a channel with dielectric boundaries and a wall normal electric field. All simulation parameters have been selected to match those in Ref. [81]. Right: comparison of electro-osmotic flow between DISCOS and the continuum electrolyte code described in Ref. [11]. The red circles and black crosses represent the DISCOS and continuum fluid velocities a short time after the electric field has been applied. The red squares and black diamonds show the flow at large time when the system has reached steady state.