Error for different <i>N</i>_<i>F</i> (<i>x</i>-axis) and <i>N</i> (<i>y</i>-axis) for simulations lasting 10⁴ time steps.

<p>The flux source was always on. The maximum difference between the simulated <i>ρ</i> and exact solution at all <i>R</i> over all 10⁴ time steps (max_<i>t</i>{<i>ɛ</i>_sim(<i>t</i>)}) is shown on a color-coded log₁₀ scale. Errors larger than 1 were truncated to 1 to make the graph more readable.</p

The functions <i>e</i>_<i>m</i>(<i>R</i>) defined in Eq (37).

<p>The hash marks above the main figure are the interpolation grid points for small <i>I</i> (Simulation 1 in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0132273#pone.0132273.t001" target="_blank">Table 1</a>).</p

Error versus simulation time for five of the six simulations listed in Table 1.

<p>The circled numbers correspond to the Simulation listed in the table. (A) The flux source is always on. (B) The flux source is randomly on for 50000 time steps. The bars on the right side are the maximum error over 50000 time steps from panel A to show that the error in the always-on simulations bounds the error of the randomly-on simulations. For all simulations, <i>N</i> = 20 and <i>N</i>_<i>F</i> = 9.</p

Parameters used in Fig 3.

<p>The circled numbers in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0132273#pone.0132273.g003" target="_blank">Fig 3</a> correspond to the simulation number (Sim.) in the Table. For each simulation, the number of near (<i>I</i>_near) and far (<i>I</i>_far) interpolation grid points are shown, as well as the sum (<i>I</i>). Simulation 6 is not shown in <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0132273#pone.0132273.g003" target="_blank">Fig 3</a> because it was the same as Simulation 5 except that it did not have the uptick in error at the end of the simulation. For all simulations, <i>N</i> = 20 and <i>N</i>_<i>F</i> = 9.</p

Algorithm for the Time-Propagation of the Radial Diffusion Equation Based on a Gaussian Quadrature

<div><p>The numerical integration of the time-dependent spherically-symmetric radial diffusion equation from a point source is considered. The flux through the source can vary in time, possibly stochastically based on the concentration produced by the source itself. Fick’s one-dimensional diffusion equation is integrated over a time interval by considering a source term and a propagation term. The source term adds new particles during the time interval, while the propagation term diffuses the concentration profile of the previous time step. The integral in the propagation term is evaluated numerically using a combination of a new diffusion-specific Gaussian quadrature and interpolation on a diffusion-specific grid. This attempts to balance accuracy with the least number of points for both integration and interpolation. The theory can also be extended to include a simple reaction-diffusion equation in the limit of high buffer concentrations. The method is unconditionally stable. In fact, not only does it converge for any time step Δ<i>t</i>, the method offers one advantage over other methods because Δ<i>t</i> can be arbitrarily large; it is solely defined by the timescale on which the flux source turns on and off.</p></div

Simulation error versus the number of nonzero elements in <i>W</i> (to represent computation time).

<p>Each point is for a different parameter set (<i>N</i>, <i>N</i>_<i>F</i>, <i>I</i>_near, and <i>I</i>_far) for a long simulation of 10⁶ time steps for the example described in the main text with the flux source always on. The arrow points to the optimal balance of high accuracy and computation speed. Errors larger than 1 were truncated to 1 to make the graph more readable.</p

<i>ρ</i>(<i>R</i>) for a randomly-on flux source at various simulation times.

<p>It illustrates both the complicated structure of these curves and to show how the interpolation grid points are more dense in the regions that need them. The points are the simulation results and the curves are the exact result. For comparison, the dashed curve is the always-on exact result at the same time step as the red curve. For clarity a very small number of interpolation grid points were used (<i>I</i> = 64).</p

Steady-State Electrodiffusion from the Nernst–Planck Equation Coupled to Local Equilibrium Monte Carlo Simulations

We propose a procedure to compute the steady-state transport of charged particles based on the Nernst–Planck (NP) equation of electrodiffusion. To close the NP equation and to establish a relation between the concentration and electrochemical potential profiles, we introduce the Local Equilibrium Monte Carlo (LEMC) method. In this method, Grand Canonical Monte Carlo simulations are performed using the electrochemical potential specified for the distinct volume elements. An iteration procedure that self-consistently solves the NP and flux continuity equations with LEMC is shown to converge quickly. This NP+LEMC technique can be used in systems with diffusion of charged or uncharged particles in complex three-dimensional geometries, including systems with low concentrations and small applied voltages that are difficult for other particle simulation techniques

Ion Correlations in Nanofluidic Channels: Effects of Ion Size, Valence, and Concentration on Voltage- and Pressure-Driven Currents

The effects of ion–ion and ion–wall correlations in nanochannels are explored, specifically how they influence voltage- and pressure-driven currents and pressure-to-voltage energy conversion. Cations of different diameters (0.15, 0.3, and 0.9 nm) and different valences (+1, +2, and +3) at concentrations ranging from 10^–6 M to 1 M are considered in 50-nm- and 100-nm-wide nanoslits with wall surface charges ranging from 0 C/m² to −0.3 C/m². These parameters are typical of nanofluidic devices. Ion correlations have significant effects on device properties over large parts of this parameter space. These effects are the result of ion layering (oscillatory concentration profiles) for large monovalent cations and charge inversion (more cations in the first layer near the wall than necessary to neutralize the surface charge) for the multivalent cations. The ions were modeled as charged, hard spheres using density functional theory of fluids, and current was computed with the Navier–Stokes equations with two different no-slip conditions

PACO: PArticle COunting Method To Enforce Concentrations in Dynamic Simulations

We present PACO, a computationally efficient method for concentration boundary conditions in nonequilibrium particle simulations. Because it requires only particle counting, its computational effort is significantly smaller than other methods. PACO enables Brownian dynamics simulations of micromolar electrolytes (3 orders of magnitude lower than previously simulated). PACO for Brownian dynamics is integrated in the BROWNIES package (www.phys.rush.edu/BROWNIES). We also introduce a molecular dynamics PACO implementation that allows for very accurate control of concentration gradients