Multiscale modeling of a rectifying bipolar nanopore: explicit-water versus implicit-water simulations

In a multiscale modeling approach, we present computer simulation results for a rectifying bipolar nanopore on two modeling levels. In an all-atom model, we use explicit water to simulate ion transport directly with the molecular dynamics technique. In a reduced model, we use implicit water and apply the Local Equilibrium Monte Carlo method together with the Nernst-Planck transport equation. This hybrid method makes the fast calculation of ion transport possible at the price of lost details. We show that the implicit-water model is an appropriate representation of the explicit-water model when we look at the system at the device (i.e., input vs.\ output) level. The two models produce qualitatively similar behavior of the electrical current for different voltages and model parameters. Looking at details of concentration and potential profiles, we find profound differences between the two models. These differences, however, do not influence the basic behavior of the model as a device because they do not influence the $z$-dependence of the concentration profiles which are the main determinants of current. These results then address an old paradox: how do reduced models, whose assumptions should break down in a nanoscale device, predict experimental data? Our simulations show that reduced models can still capture the overall device physics correctly, even though they get some important aspects of the molecular-scale physics quite wrong; reduced models work because they include the physics that is necessary from the point of view of device function. Therefore, reduced models can suffice for general device understanding and device design, but more detailed models might be needed for molecular level understanding.


Introduction
Our devices are getting smaller and smaller. Nanodevices that involve ion transport generally have regions that can be considered macroscopic from the molecular point of view, with characteristic size of µm scale and larger (e.g., wide part of a conical nanopore, access regions, bulk regions with electrodes). Continuum theories are often sufficient to model and compute devices built of such regions. These reduced (or coarse-grained) models make approximations like pointcharge ions and implicit water as a trade-off for fast computation time. As the dimensions of the devices shrink, however, nanoscopic regions appear (e.g., narrow bottlenecks of nanopores), where the assumptions of these reduced models are no longer valid and molecular models and particle simulations like all-atom molecular dynamics (MD) are needed to reveal molecular mechanisms. To connect with experimental data then requires a multiscale modeling approach, where different models/methods with different resolutions are applied to different regions. In our approach, multiscaling is applied between models for the whole system in separate simulations rather than between separate regions in a single model and simulation. The multiscaling part of our procedure is connecting the various modeling levels; molecular insights obtained from atomic models, for example, can be used to design better reduced models.
Unexpectedly, however, reduced models often reproduceand predict-experimental data, even though they should, in principle, not work because their underlying atomic physics is too approximate; one should not assume that water molecules are absent or that ions are point charges when the pore they go through is only a few times wider than the "real" ion or water. Yet, the Poisson-Nernst-Planck (PNP) theory can reproduce experimental data for rectifying nanopores [1,2,3,4,5,6]. Other examples include reduced models using hard sphere ions reproducing and predicting nonlinear phenomena in biological ion channels [7,8,9,10,11,12,13,14] and in nanopores [15,16]. Moreover, in many cases including nanopores and ion channels, reduced models are the only way to connect with experimental results, as the required low ion concentrations and small applied voltages are inaccessible to all-atom MD simulations.
The goal of this work is to report computational results for a nanodevice, a rectifying nanopore, on two distinct modeling levels that differ in the treatment of water. In one model, water is a continuum background and the ions are charged, hard spheres that move via drift-diffusion. In the other model, MD simulations with explicit water molecules are used. This simple case study allows us to identify the effects of using an implicit water model versus having corpuscular water molecules.
Consistent with the work cited above, we find that both models produce qualitatively the same input-output relations (e.g., current vs. voltage, pore surface charge, or ion concentration). On the other hand, the molecular-scale picture inside the nanopore for the two models is quite different, as one would expect. The solution to this paradox is that both models produce qualitatively the same axial ion concentration profiles (i.e., cross-section averaged profiles), and because these are the first-order determinants of current, the two models produce the same device characteristics.
Our results suggest that reduced models that reproduce experimental results are, in fact, valid because they capture the overall device physics correctly despite the incomplete description of the molecular level phenomena. Therefore, they are useful for understanding the device-level physics that produces the input-output relations and for device design. At that level, the molecular details are of second order importance. However, to understand the physics at an atomic scale (e.g., ion and water structure inside a confining nanopore), more detailed molecular models must be used.

Models
In an all-atom model, water molecules are modeled explicitly and the system is simulated with the MD technique. The ionic current, driven by an external electric field, is simulated directly by counting the ions that passed through the pore (more detail is found in the Appendix and in the ESI).
Coarse-grained models reduce the number of degrees of freedom in the Hamiltonian by integrating out some of these degrees of freedom and replacing them with a response function, or simplify model potentials describing certain parts of the system. In this work, the reduced model describes water as an implicit continuum background, where the effect of water is replaced by various response functions: (1) The ability of water to affect the movement of ions through friction (a dynamic effect) is taken into account by a diffusion coefficient, D i (r), in the Nernst-Planck (NP) transport equation with which we compute the flux: where j i (r) is the particle flux density of ionic species i, c i (r) is concentration, k is Boltzmann's constant, T is temperature, and µ i (r) is the electrochemical potential profile. Its gradient is the driving force of the steady-state transport.
(2) The ability of water to screen the charges of ions (an energetic effect) is taken into account by a dielectric constant, ε, in the denominator of the Coulomb-potential acting between the charged hard spheres with which we model the ions: where q i is the charge and R i is the radius of ionic species i, ε 0 is the permittivity of vacuum, and r is the distance between the ions. We simulate this model with the Local Equilibrium Monte Carlo (LEMC) technique [17,18,19], which is practically a grand canonical Monte Carlo simulation devised for a nonequilibrium situation. The input variable of the LEMC simulation is the chemical potential profile, µ i (r), that is not constant for a system out of equilibrium, but a space-dependent quantity. The output variable is the concentration profile, c i (r). Thereby, the LEMC simulation establishes the relation between c i (r) and µ i (r) necessary to apply the NP equation. The LEMC method correctly computes volume exclusion and electrostatic correlations between ions, so it is beyond the mean-field level of the PNP theory that is routinely applied for nanopores.
The NP equation is solved iteratively with the LEMC simulations in an iteration procedure that ensures that the continuity equation, ∇ · j i (r) = 0, is satisfied (more detail is found in the Appendix and in the ESI). The resulting NP+LEMC technique provides a solution for the statistical mechanical problem (e.g., the c i − µ i relation) on the basis of particle simulations, while it still gives an approximate indirect solution for the dynamical problem through the NP equation. Direct simulation of ionic transport in the implicit-water framework is commonly done by the Brownian Dynamics method [20,21,22]. The advantage of the NP+LEMC technique is that we can easily handle cases where direct sampling of ions crossing the pore is problematic such as the case of the bipolar nanopore, through which currents are small due to the depletion of ions in the oppositely charged zones.

Questions
The questions that we pose in this study are the following. Is the implicit-water model a good approximation to the explicit-water case from the device point of view? Do the two models provide the same device behavior, namely, have the same transfer function for the device? The transfer function describes the relation between the control signal (input) and the response (output) that the model gives to the input signal. Strictly speaking, this is the electrical current as a function of voltage that drives the current. We can define transfer functions, however, in a more general way: the electrical current as a function of experimentally tunable parameters such as pore charge and bulk concentrations.
The transfer functions are the same on the two modeling levels because we study the same device with the two models. The implicit model is justified on the basis of its ability to reproduce the transfer functions given by the detailed model. The MD model, therefore, is the gold standard in this comparison because it contains all the molecular details (explicit water) that are missing from the implicit-water model. If we consider the device as a black box model, we are satisfied with just answering this question. As the dimensions of the devices are getting smaller, however, we want to open the box and look inside to understand how it works on the level of molecules, which is necessary to design new devices.
Even if the explicit-and implicit-water models work similarly on the device-level, opening the box reveals that there are profound differences between the two models due to the presence or absence of explicit water molecules. In this case, new questions arise. How can the transfer functions be similar despite the differences? Specifically, why can water molecules be averaged into an implicit background even in confined spaces? Using a reduced model we may be ignoring crucial details. Looking at too much detail, however, can result in missing the forest for the trees. Unknown errors from inadequately calibrated parameters of detailed models may produce ambiguous artifacts. What are the details that cannot be coarse-grained and should be included in the reduced model to reproduce the device function? Bipolar nanopore models To address these questions, we have chosen a bipolar nanopore as a test system because it has a very characteristic transfer function: the pore rectifies electrical current. The output is the current given as a response to the input, the voltage. The current is markedly different for different signs of the voltage due to the asymmetrical surface charge distribution on the pore wall. The sign of the surface charge changes from positive to negative along the central axis of the pore in this study -we call these N and P regions, respectively, according to the nomenclature of the semiconductor literature ( Fig. 1). In the N zone (with positive surface charge), the anions are the counter-ions and therefore the majority charge carriers, while the cations are the co-ions and the minority charge carriers. In the P zone, the situation is reversed. Pore regions with opposite surface charges can be achieved by chemical modifications in the case of PET nanopores by transforming carboxyl groups into amino groups by a coupling agent [23], for example. This modification can even be reversible by making the chemical moieties redoxsensitive [24]. The pore charges in the N and P zones do not drive ion tranport; rather, they modulate the concentrations of the various ionic species in the respective zones. The driving force of drift-diffusion of ions is the difference of the concentrations and/or applied voltage on the two sides of the membrane [25,26,27,28,29,30,31,32,33,34,1,35,36,37,15,16,38,39].
Fabrication of such synthetic nanopores have become possible recently [40,41,42,43,44,45,46,47]. Their common property (that distinguishes them from micropores) is that their radii are smaller than (or similar in size to) the characteristic screening length of the electrolyte (NaCl in this study) that fills the channel. The dimensionality of these systems makes them central building blocks of devices for various applications from DNA sequencing [48] through switches [49] to chemical and biosensing [50,51,52,53]. Rectification appears (and the nanopore can be called a nanofluidic diode) if the nanopore is asymmetrical in geometry (e.g. conical nanopores) or in charge distribution (e.g. bipolar nanopores studied in this work). What makes the bipolar nanopore more suitable to be the test system of our multiscale modeling study is the fact that the mechanism of rectification is mainly electrostatic in nature and that it rectifies even if the dimensions of the pore are small (1 nm radius and 6 nm length in this work) and therefore suitable to study with MD.
In this study we pay special care to modeling the nanopore in the MD and the NP+LEMC calculations as similarly as possible, so that the major difference between the two systems is the treatment of water, compared to which other dissimilarities are minor. This makes this work a case study for studying the differences between the explicit and implicit water models. We connect the two modeling levels by estimating the diffusion coefficients of ions in the pore (for the NP+LEMC calculations) from the MD simulations. We fit the parameters of the D i (r) functions so that NP+LEMC results for the ionic currents reproduce MD results as closely as possible. In this study, we choose the simplest assumption for the D i (r) profile: it is a piecewise constant function that is different inside the pore and in the bulk: D values accounts for all those effects that are different in the pore and in the bulk including different structure of water molecules around the ions and the presence of the pore (confinement and being trapped in electrostatic wells created by the pore charges). We fit the D pore i values for one state, fix them, and check whether they are appropriate in other cases too.
We work with a cylindrical nanopore model whose dimensions are relatively small, only ∼ 1 nm in radius and ∼ 6 nm in length. This radius is quite close to that attainable in experiments at the tip of conical nanopores. The length is far from the experimentally realistic µm length scale, but we are restricted by the fact that we use particle simulations. The small pore radius allows the overlapping of the double layers formed at the pore walls at the 0.5 and 1 M concentrations used in this study (wider pores are considered elsewhere [66]). Large concentrations were necessary to achieve a reasonable sampling in the MD simulations. Small pore radius is also useful for the purpose of this study because it amplifies the role of water molecules.
The surfaces of the pore and the membrane are smooth. This is advantageous because we can evaluate the results more easily, we can avoid artifacts from surface roughness, and we can create toy models for both the MD simulations (using a carbon nanotube (CNT) as template, see Illustration of the geometry and the rectification mechanism. Top panels: Nanopore structure made by the VMD package [73]. The pore is formed by a carbon nanotube (CNT) between two carbon nanosheets (CNS) defining the membrane. Distance of CNS atoms in the two sheets (width of membrane) is H = 6.035 nm, while the distance of CNT atoms from the pore axis (pore radius) is R c = 1.136 nm. The pore is made charged by placing partial charges on the C atoms of the CNT. The charges are positive on the left hand side (blue dots, N region) and negative on the right hand side (red dots, P region) of the pore. The value of partial charges is chosen so that the average surface charge density corresponds to a prescribed value, ±σ . The series of spheres inside the pore represent trajectories of Cl − ions through the pore plotted in 10 ps time intervals. The Cl − ion was randomly chosen of those that crossed the pore. The trajectories are from simulations for c = 1 M NaCl and σ = 1 e/nm 2 at ±200 mV voltages (200 mV is ON -left panels, while -200 mV is OFF -right panels). Middle panels: randomly chosen snapshots (from a video available in the ESI) with blue and red spheres representing Na + and Cl − ions, respectively. The depletion zone of Na + is formed in the entire N zone, not only in the PN junction in the middle (the same is true for Cl − in the P region). The depletion zones are deeper in the OFF state. Bottom panels: illustrative MD concentration profiles (more detailed profiles are shown in Fig. 3).
NP+LEMC calculations (using hard walls, see Fig. SI 12). This way, the results of the two methods are comparable and the main difference between them is the treatment of water.

Results and Discussion
Results for the device level First, we looked at the nanopore from the device perspective, namely, we studied how the response function behaves in the two modeling frameworks. Because MD sampling is weak for small concentrations and voltages, we restrict ourselves to voltages ±200 mV and concentration c = 1 M in the main text. The effect of the bulk concentration and voltage-dependence (current-voltage curves) is found in the ESI showing mainly NP+LEMC results with a few MD data (Figs. SI 1 and SI 5). For a full comparison between the two models, we focus on the dependence of currents on the surface charge density on the nanopore wall, σ . We place σ surface charge on the left hand side of the pore wall and −σ on the right hand side (see Fig. 1), therefore, the value of σ characterizes the strength of the polarity of the pore. We define rectification as the ratio of the magnitudes of the ON (200 mV) and the OFF (-200 mV) currents (ground is on the left hand side). values transferable? If they are transferable, the reduced model "knows" the basic physics necessary to describe the device behavior, so no additional fitting is needed. This means that the approximations built into the reduced model smear only the "unimportant" degrees of freedom into response functions (ε and D i (r)), but leave the "important" degrees of freedom unaffected. A crucial element of multiscaling is identifying the "important" degrees of freedom by doing the calculations at the various resolutions and comparing the results.
Let us examine Fig. 2 in more detail. When the pore is uncharged (σ = 0 e/nm 2 ), the currents at the positive and negative voltages are the same (apart from statistical uncertainties; see ESI for details), because the pore is symmetric. When σ is increased, the currents in the ON and OFF states increasingly deviate. The current (in absolute value) increases in the ON If we want to look into the black box, we need to analyze space-dependent concentration and potential profiles. From the point of view of current, which is the integral of the flux density over the cross section of the pore, the z-dependent concentration profiles (averages over cross sections) are the most relevant profiles. They describe which ion is present in which region (N and P) of the nanopore and are the first-order determinants of the current.
The top panels of Fig. 3 show concentration profiles obtained from MD (left panels) and NP+LEMC (right panels) calculations. The two models describe the mechanism of recti- In the bottom panels, the ratio of concentration profiles in the ON and OFF states are shown for various values of σ (see the numbers near the curves). Concentration profiles have been computed in the same way in the two models: the average number of ions in a slab has been divided by the effective volume available for the ions. Inside the pore, the same cross section was used (radius 1 nm) in both cases. fication in qualitative agreement. As a primary effect, the pore charges produce depletion zones for the respective ions, e.g., the positive pore charges on the left hand side (N zone) result in depletion zones in the Na + profiles (the same is true for the Cl − profiles in the P zone). The Na + ions, therefore are co-ions in the N zones, while the Cl − ions are counter-ions there (reversed in the P zone). In this system, diffusion is limited for both ionic species by their respective depletion zones because these zones are the largest resistance elements if we imagine the ionic pathway as resistors connected in series in an equivalent circuit; if an ion species is not there, it cannot conduct. Depletion zones, therefore, are more important than peaks.
As a secondary effect, the applied field modulates the effect of pore charges. It increases the concentrations (of both ions) in the ON state, while it decreases them in the OFF state. Depletion zones of the ions, therefore, are deeper in the OFF state than in the ON state. This is the basic reason of rectification. It is important to stress, however, that the peaks are also lower in the OFF state. The depletion zones are not independent from the peaks. Co-ions and counter-ions are strongly correlated electrostatically, therefore, co-ions are present in their depletion zones because the counter-ions are present there.
As far as the σ -dependence is concerned, ON concen-trations increase, while OFF concentrations decrease with increasing σ as shown by the bottom panels of Fig. 3, where the ratio of the ON and OFF concentrations is shown for different σ values (the concentration profiles themselves are in Fig. SI 4). The ratio increases with increasing σ . The ratio is larger in the case of the MD data (corresponding to deeper depletion zones), which is reflected in the smaller OFF currents compared to the NP+LEMC results. Better agreement between MD and NP+LEMC could be achieved by adjusting different diffusion constant values in the N and P zones (different mobilities in the N and P zones are implied by the top panels of Fig. 1), but this is beyond the scope of this study. The results obtained for the operation of the nanopore as a device (and the accompanying c i (z) profiles) are quite similar on the two modeling levels. This is surprising given that the treatments of water in the two models are so profoundly different. Consequences of the differences can be revealed if we open the black box even wider and dive into details beyond the cross-sectionally averaged c i (z) profiles.

Results for the molecular level
The presence of water molecules and their effect on ions can be seen by plotting the full c i (z, r) profiles. Figure 4A shows contour plots for the ON state as obtained from MD simulations (the plot for the OFF state is seen in Fig. SI 7). A clear layering structure as a function of the r coordinate is present. There are two distinct peaks of water, one is near the pore wall solvating the pore charges. The ions' behavior depends on the region in which we observe them. Cl − ions have a large peak just "behind" the solvating water layer in the N region, while they rather accumulate in the pore center in the P region. Na + ions, on the other hand, show the opposite behavior, exhibiting two peaks in the P region. These effects are better (and more quantitatively) observed in Fig. 4B, where the c i (r) profiles are plotted averaged over the z-coordinate in the N region (top panel) and in the P region (bottom panel). Here, the NP+LEMC profiles for Na + and Cl − are also shown. Although they are unstructured, they reproduce the relative quantities of the two ion species in the respective regions.
Therefore, despite the differences in the radial structure, the MD and NP+LEMC ionic profiles show basically the same behavior in the z dimension. This is the behavior that is relevant for the calculation of the current, and therefore, for the function of the device. The r-dependence seems to belong to the class of the "unimportant" degrees of freedom vis-à-vis ion current, because its accurate reproduction is not necessary to properly model the operation of the nanopore as a rectifying diode. The z-dependence, on the other hand, is crucially important. This explains why 1-dimensional PNP calculations work so well for nanopores. The water molecules then also seem to belong to the class of "unimportant" degrees of free-dom, in the sense that they do not contribute to the current (they have no net charge) and do not affect the axial c i (z) profiles of ions substantially. Also, the treatment of their dynamic effect on the ions as a friction via the diffusion coefficients in the NP equation seems to be a sufficient approximation.
Their other effect on ions, screening, also works very differently in the implicit and explicit-water models. To perceive these profoundly different ways of screening as given by MD and NP+LEMC simulations, we consider the electrical potential profiles produced by the various species in the two models. The electrostatic potential from water is present in the MD simulations explicitly and provides screening of the potential from ions (and pore charges) in an additive way. In the NP+LEMC simulations, screening is provided by polarization charges (1/ε − 1)q that are right on top of the q ionic point charges. In effect, this correponds to just damping the electric field of ions by dividing it by ε (see Eq. 2). Figure 5 shows potential profiles as contour plots over the (z, r) coordinates obtained from MD (left panel) and NP+LEMC (right panel) at the same conditions. These potentials do not contain the applied potential, only those produced by free charges in the system: ions, pore charges, and water (if present). Top halves are for the ON state, while bottom halves are for the OFF state. The differences between the two models are striking.
In the implicit-water model (Fig. 5B), only ions can respond to changes in the electric field. The (1/ε − 1)q polar- ization charges move together with the ions so they cannot respond to the electric field independently. The effect of the applied field, therefore, is clearly visible. More ions are available for screening in the ON state, so the potential is smaller (in absolute value) in the ON state. In the OFF state, due to the low concentration of ions, the pore charges dominate the potential.
In the case of the explicit-water model (Fig. 5A), on the other hand, water molecules also produce a counter-field to external effects. The first-order external effects to which both ions and water respond are the applied field and pore charges (these are fixed). Water molecules also respond to the field of ions (and vice versa, of course). Because ions move slowly and they alone (without water) form a low-density (dilute) system, the characteristic speed at which they form the ionic distribution is relatively small. Water molecules, on the other hand, are high-density (liquid) and rotate quickly, so they accomodate themselves to the movement of ions and other external fields fast. The electrolyte, therefore, is a conducting material in which mobile charges are found in high density that respond to external fields and exert a counter-field that produce close to constant potential in the electrolyte.
Because the potentials in the implicit-water model are uniformly damped by ε = 78.5, the potential varies in a narrow range (−2.5 < eΦ(z, r)/kT < 2.5). In the explicit-water case, on the other hand, the total potential changes in a much wider range (Fig. 5B vs. A). For this reason, the effect of the relatively small applied field (±200 mV) is hardly visible.
The effect of the applied field, on the other hand, is better visible on the potentials produced by the ions (together with pore charges) and water molecules separately. These poten-tials are plotted in Fig. 6. Figures 6A and B show the results for the ON and OFF states, respectively. Water molecules produce an electrical potential (bottom halves) that has an opposite sign compared to the potential of ions and pore charges (top halves).
Differences in the potentials of ions and pore charges in the ON and OFF states are consequences of the differences in the ionic profiles in the ON and OFF states. Water molecules respond easily to the changes in ion distributions (as described above), so their potentials are also different in the ON and OFF states.
Another important difference between the explicit-and implicit-water models is whether ionic double layers are formed near the membrane on the two sides of it. In the implicit-water case they are shown by the splitting cation and anion profiles (Fig. 3 and Fig. SI 9) and by the gradually rising/declining potential profiles ( Fig. 5B and Fig. SI 10) on the two sides of the membrane. The signs of the double layer potentials are the opposite in the ON and OFF states and their formation is closely correlated with the relative quantities of counter-ions vs. co-ions in the pore region adjacent to the double layer (see our other papers [65,66] for details). The presence of the double layers, however, is not necessary to produce rectification. They are absent in the explicit-water MD simulations, still, the different currents in the ON and OFF states are present.

Conclusions
We found definite differences between the results of the explicitand implicit-water models with respect to (1) how they de- scribe the radial dependence of the particle profiles, (2) whether the double layers on the two sides of the membrane are formed or not, and (3) how the screening of ionic and pore charges by water is done. Still, the z-dependence of the ionic concentration profiles is similar in the two models. The effect of the pore charges producing the depletion zones and the effect of the applied field in modulating the effect of pore charges are so powerful that ions respond to them in a similar way no matter whether the water molecules are there explicitly or their screening effect is just a scaling factor in the denominator of the Coulomb potential.
While this is a simple case study, it does answer the paradox of why reduced models work (i.e., reproduce experimental results) in regimes (e.g., confined geometries) where they ought not to work because their molecular-level physics is too approximate. Our results show that some atomic details do not matter for device input-output properties, and that reduced models that reproduce experiments get the overall device physics right, even though they may get some aspects of the molecular-scale physics wrong. That is, reduced models show what device-level physics is necessary/sufficient to predict device characteristics and to design new devices. The results of reduced models should not be overanalyzed, however. The nanoscale-level results they produce are likely to be incorrect in many ways and MD simulations should be used if these details are desired.
For the foreseeable future, however, fully atomistic simulations will probably remain too computationally expensive to regularly predict experiments, and reduced models will remain the only way to compare with experiments in many cases. Our work provides computational insight into how these reduced models should be interpreted.

Appendix
All-atom MD simulations were performed with the GRO-MACS (v.5.0.4) program suite [74,75] using the leap-frog integrator with a 2 fs time step. The CHARMM27 force field [76] was implemented that included ∼11,000 SPC water molecules, 110-190 Lennard-Jones (LJ) and point charge type anions/cations (depending on the concentration and pore charge density), and LJ atoms for the CNT and CNS. The dimensions of the simulation box were 6×5.2×16.8 nm, with periodic boundary conditions (PBC) applied in all spatial directions. The systems were thermostated to 298.15 K by a modified version of the Berendsen (velocity rescaling) algorithm [77]. An ion was considered to cross the channel if it is initially at one side of the membrane and then ends at the opposite side of the membrane after propagating through the channel.
In the reduced model designed for NP+LEMC calculations the membrane and the pore penetrating it are defined by hard walls with which the ions cannot overlap. The cylindrical pore's radius and length were calculated to mimic the CNT model of the MD simulations as closely as possible on the basis of an estimated distance of closest approach of ions to the carbon atoms. We used the values R = 0.97 nm and H = 6.4 nm for the pore radius and length, respectively. The fractional point charges have been placed at the same positions as in the CNT model. The Na + and Cl − ions are charged hard spheres with radii 0.095 and 0.181 nm) immersed in a dielectric continuum (ε = 78.5) that models the solvent implicitly (Eq. 2). The bulk diffusion coefficients of Na + and Cl − were D bulk Na + = 1.333×10 −9 m 2 s −1 and D bulk Cl − = 2.032×10 −9 m 2 s −1 , respectively. The problem is solved on a discretized grid iteratively; the electrochemical potential is changed until the continuity equation, ∇ · j i (r) = 0, is satisfied. The resulting profiles fluctuate around limiting distributions due to statistical uncertainties in the LEMC simulations. The final results are obtained as running averages. Computational time is measured in days (sometimes, weeks) for MD simulations, in hours for the NP+LEMC simulations, but only in seconds for PNP calculations.
Boundary conditions are treated differently in MD and NP+LEMC. In MD, a ±0.012 V/nm homogeneous external electric field was applied along the z-direction to achieve ±200 mV potential difference. Ions leaving the cell at one end are fed back on the other end due to PBC applied in the z dimension. In NP+LEMC, the simulation cell is a finite cylinder due to the rotational symmetry [17,18,19]. Two cylindrical compartments on the two sides of the membrane represent the two bulk regions between which the ion transport flows. Both concentrations and electrical potentials are prescribed on the half-cylinders confining these bulk regions.
The potential profiles can be computed either (1) by inserting test charges and sampling the potential on the fly, or (2) by solving Poisson's equation for the averaged (z, r)-dependent charge profiles obtained from the simulations (with Dirichlet boundary conditions of assuming zero potential on the closed surface of the confining cylinder). Agremeent between the two methods was excellent for the NP+LEMC data, while it was reasonable for the MD data (see Fig. SI 8). The potential profiles shown in Figs. 5 and 6 have been obtained from the second method.