^{1}, R. E. Jones

^{2}, B. J. Debusschere

^{2}and O. M. Knio

^{3,a)}

### Abstract

In this article, uncertainty quantification is applied to molecular dynamics (MD) simulations of concentration driven ionic flow through a silica nanopore. We consider a silica pore model connecting two reservoirs containing a solution of sodium (Na +) and chloride (Cl−) ions in water. An ad hoc concentration control algorithm is developed to simulate a concentration driven counter flow of ions through the pore, with the ionic flux being the main observable extracted from the MD system. We explore the sensitivity of the system to two physical parameters of the pore, namely, the pore diameter and the gating charge. First we conduct a quantitative analysis of the impact of the pore diameter on the ionic flux, and interpret the results in terms of the interplay between size effects and ion mobility. Second, we analyze the effect of gating charge by treating the charge density over the pore surface as an uncertain parameter in a forward propagation study. Polynomial chaos expansions and Bayesian inference are exploited to isolate the effect of intrinsic noise and quantify the impact of parametric uncertainty on the MD predictions. We highlight the challenges arising from the heterogeneous nature of the system, given the several components involved, and from the substantial effect of the intrinsic thermal noise.

This work was supported by the U.S. Department of Energy Office of Science through the Applied Mathematics program in the Office of Advanced Scientific Computing Research (ASCR) and the Laboratory Directed Research and Development (LDRD) program at Sandia National Laboratories. Sandia is a multiprogram laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Corporation, for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-AC04-94AL85000.

We would also like to acknowledge helpful discussions with P. Crozier (Sandia) and H. Adalsteinsson (currently at Google) regarding the concentration control algorithm.

I. INTRODUCTION

II. ATOMISTIC SYSTEM AND MD SIMULATIONS

A. Atomistic system

B. Force-field

C. MD simulations

1. Equilibration

2. Concentration control

III. DEPENDENCE ON THE PORE DIAMETER

A. Velocity profiles

B. Ionic conductance

IV. SENSITIVITY TO THE GATING CHARGE

A. UQ formulation

B. Bayesian regression

1. Collection of observations

2. Noise model and likelihood function

C. Bayesian regression: Results

1. Bayes factor

2. Posterior predictive check and predictive uncertainty

V. CONCLUSIONS

### Key Topics

- Sodium
- 74.0
- Silica
- 29.0
- Nanoporous materials
- 21.0
- Ionic conduction
- 16.0
- Molecular dynamics
- 13.0

##### B82B1/00

## Figures

(a) Schematic of the dual-reservoir nanopore system. The silica is a rectangular block with a cylindrical pore of nominal radius R. Only few surface layers of the silica block have dynamics (gray), shown in panel (b), while the rest comprises the “frozen” atoms (green). The volume of fluid is nominally .

(a) Schematic of the dual-reservoir nanopore system. The silica is a rectangular block with a cylindrical pore of nominal radius R. Only few surface layers of the silica block have dynamics (gray), shown in panel (b), while the rest comprises the “frozen” atoms (green). The volume of fluid is nominally .

Silica pore model of dimensions (x, y, z) L x ≈ 54, L y ≈ 60, and Å obtained from a bulk crystal silica structure after carving a pore of nominal radius R, ensuring stoichiometry, and hydroxylating the surface and pore defects. Color legend: silicon Si (ochre), bulk oxygens O bulk (blue), hydroxide oxygens O hy (yellow), and hydroxide hydrogens H hy (black).

Silica pore model of dimensions (x, y, z) L x ≈ 54, L y ≈ 60, and Å obtained from a bulk crystal silica structure after carving a pore of nominal radius R, ensuring stoichiometry, and hydroxylating the surface and pore defects. Color legend: silicon Si (ochre), bulk oxygens O bulk (blue), hydroxide oxygens O hy (yellow), and hydroxide hydrogens H hy (black).

Snapshots of cross-sectional views of the systems for D = 12.5 (a), 17 (b), 21 (c), and 27 Å (d), taken during the steady state of the CC stage. Color legend: Si (ochre), O bulk (blue), O hy (yellow), H hy (black), O w (red), H w (white), Na+ (purple), Cl− (cyan). The silica is nearly transparent for visualization convenience.

Snapshots of cross-sectional views of the systems for D = 12.5 (a), 17 (b), 21 (c), and 27 Å (d), taken during the steady state of the CC stage. Color legend: Si (ochre), O bulk (blue), O hy (yellow), H hy (black), O w (red), H w (white), Na+ (purple), Cl− (cyan). The silica is nearly transparent for visualization convenience.

(a) Time/bin-averaged radial profile of the axial velocity v z , for Na+ and Cl− computed for one replica of the case D = 21 Å during the steady phase of the CC stage. (b) Time/bin-averaged axial velocity v z plotted as a function of the z-coordinate computed for one replica of D = 21 Å at the final step of the CC stage for water, Na+ and Cl−. In (a), for the sake of visualization, the data obtained from the spatial binning are interpolated over a finer mesh. In (b) the filled circles are plotted at the location of the 24 bins.

(a) Time/bin-averaged radial profile of the axial velocity v z , for Na+ and Cl− computed for one replica of the case D = 21 Å during the steady phase of the CC stage. (b) Time/bin-averaged axial velocity v z plotted as a function of the z-coordinate computed for one replica of D = 21 Å at the final step of the CC stage for water, Na+ and Cl−. In (a), for the sake of visualization, the data obtained from the spatial binning are interpolated over a finer mesh. In (b) the filled circles are plotted at the location of the 24 bins.

Time evolution of the running average conductance, G(t), for Na+ (a) and Cl− (b) plotted for all 5 replicas and each diameter value showing the variation in steady state values and the timescale at which steady values are achieved.

Time evolution of the running average conductance, G(t), for Na+ (a) and Cl− (b) plotted for all 5 replicas and each diameter value showing the variation in steady state values and the timescale at which steady values are achieved.

(a) Steady-state replica values (markers) of the Na+ (blue) and Cl− (red) conductance as a function of the nominal pore diameter D, with the superimposed lines outlining the replica-averaged values. (b) Coefficient of variation, σ/μ, i.e., the ratio of the standard deviation, σ, over the mean, μ, of the conductance data plotted as function of the pore diameter, D, for both ions.

(a) Steady-state replica values (markers) of the Na+ (blue) and Cl− (red) conductance as a function of the nominal pore diameter D, with the superimposed lines outlining the replica-averaged values. (b) Coefficient of variation, σ/μ, i.e., the ratio of the standard deviation, σ, over the mean, μ, of the conductance data plotted as function of the pore diameter, D, for both ions.

Snapshots obtained for D = 12.5 (a), 17 (b), 21 (c), and 27 Å (d), during the steady state of the CC stage, showing the distribution of water molecules around the ions passing through the pore: Na+ is color-coded purple, while Cl− is color-coded cyan.

Snapshots obtained for D = 12.5 (a), 17 (b), 21 (c), and 27 Å (d), during the steady state of the CC stage, showing the distribution of water molecules around the ions passing through the pore: Na+ is color-coded purple, while Cl− is color-coded cyan.

Panel (a) Data set of the ionic conductance, G, computed for all replicas as a function of the gating charge density q surf for Na+ and Cl−. The superimposed solid lines connect the replica-averaged values. Panel (b) The corresponding variances plotted as function of q surf with a logarithmic scale on the y-axis. Note that since q surf = 0.265619083 C/m2 yields vanishing Na+ conductance for all three replicas, the corresponding variance is not shown in panel (b).

Panel (a) Data set of the ionic conductance, G, computed for all replicas as a function of the gating charge density q surf for Na+ and Cl−. The superimposed solid lines connect the replica-averaged values. Panel (b) The corresponding variances plotted as function of q surf with a logarithmic scale on the y-axis. Note that since q surf = 0.265619083 C/m2 yields vanishing Na+ conductance for all three replicas, the corresponding variance is not shown in panel (b).

(Top row) Scatter plots of 20 000 MCMC chain samples obtained for π(g 0, g 1) for Na+ (a) and Cl− (b), extracted from the original MCMC chain by removing the burn-in period comprising the first 15 000 samples. (Bottom row) Marginalized posterior, π(g k ), of each PC coefficient g k , k = 0, …, 4, obtained via KDE for Na+ (c) and Cl− (d). All panels show the results for a linear (P = 1), quadratic (P = 2), cubic (P = 3), and quartic (P = 4) expansion.

(Top row) Scatter plots of 20 000 MCMC chain samples obtained for π(g 0, g 1) for Na+ (a) and Cl− (b), extracted from the original MCMC chain by removing the burn-in period comprising the first 15 000 samples. (Bottom row) Marginalized posterior, π(g k ), of each PC coefficient g k , k = 0, …, 4, obtained via KDE for Na+ (c) and Cl− (d). All panels show the results for a linear (P = 1), quadratic (P = 2), cubic (P = 3), and quartic (P = 4) expansion.

Panel (a) shows the comparison between the data points of the Na+ (black circles) conductance and the corresponding predictions obtained from the MAP estimate of the PC regression function, M(ξ), computed for a linear (P = 1), quadratic (P = 2), cubic (P = 3), and quartic (P = 4) expansion as a function of q surf . The corresponding results for Cl− are shown in panel (b).

Panel (a) shows the comparison between the data points of the Na+ (black circles) conductance and the corresponding predictions obtained from the MAP estimate of the PC regression function, M(ξ), computed for a linear (P = 1), quadratic (P = 2), cubic (P = 3), and quartic (P = 4) expansion as a function of q surf . The corresponding results for Cl− are shown in panel (b).

Panel (a) shows the data-based variance (black circles) of the Na+ conductance versus the corresponding predictions obtained from the MAP estimate of the noise PC coefficients {d 0, d 1}, as a function of the order, P, of the regression function, M(ξ). The corresponding results for Cl− are shown in panel (b). Both plots are presented with a log scale on the y-axis.

Panel (a) shows the data-based variance (black circles) of the Na+ conductance versus the corresponding predictions obtained from the MAP estimate of the noise PC coefficients {d 0, d 1}, as a function of the order, P, of the regression function, M(ξ). The corresponding results for Cl− are shown in panel (b). Both plots are presented with a log scale on the y-axis.

Results showing the posterior predictive check samples (gray) obtained for the Na+ (a) and Cl− (b) conductance using a third-order and fourth-order PC representation, respectively. Superimposed to the plots, we report the original data color-coded blue for Na+ and red for Cl−, and the mean (black square) and error bars for , where is the standard deviation calculated from the original data set of conductances used in the inference.

Results showing the posterior predictive check samples (gray) obtained for the Na+ (a) and Cl− (b) conductance using a third-order and fourth-order PC representation, respectively. Superimposed to the plots, we report the original data color-coded blue for Na+ and red for Cl−, and the mean (black square) and error bars for , where is the standard deviation calculated from the original data set of conductances used in the inference.

Results of the “model uncertainty” analysis obtained for Na+ (a) and Cl− (b) showing how the posterior uncertainty in the inferred PC representations of the conductances is reflected in the corresponding predictions.

Results of the “model uncertainty” analysis obtained for Na+ (a) and Cl− (b) showing how the posterior uncertainty in the inferred PC representations of the conductances is reflected in the corresponding predictions.

## Tables

Lennard-Jones parameters (ɛ, σ) and Coulombic charge, in multiples of the electron charge |e|, for each atom type present in the system: {H w , O w } are water hydrogen and oxygen atoms, respectively, {H hy , O hy } are hydrogen and oxygen atoms appearing in an hydroxide group, O bulk are oxygen atoms in the bulk of the silica, Si are silicon atoms, and {Na+, Cl−} are the salt ions. The Lorentz-Berthelot mixing rules , and are used to define the interspecies Lennard-Jones interactions. For the remaining set of force-field parameters, refer to Ref. 38 for silica, to Refs. 42,44 for water, and to Ref. 45 for the ions.

Lennard-Jones parameters (ɛ, σ) and Coulombic charge, in multiples of the electron charge |e|, for each atom type present in the system: {H w , O w } are water hydrogen and oxygen atoms, respectively, {H hy , O hy } are hydrogen and oxygen atoms appearing in an hydroxide group, O bulk are oxygen atoms in the bulk of the silica, Si are silicon atoms, and {Na+, Cl−} are the salt ions. The Lorentz-Berthelot mixing rules , and are used to define the interspecies Lennard-Jones interactions. For the remaining set of force-field parameters, refer to Ref. 38 for silica, to Refs. 42,44 for water, and to Ref. 45 for the ions.

Computed values of , obtained for all four different models and each ion Na+ and Cl−. Note that due to the logarithmic scale adopted, for a given ion, the corresponding matrix of values is antisymmetric.

Computed values of , obtained for all four different models and each ion Na+ and Cl−. Note that due to the logarithmic scale adopted, for a given ion, the corresponding matrix of values is antisymmetric.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content