Uncertainty quantification in MD simulations of concentration driven ionic flow through a silica nanopore. Part II: uncertain potential parameters. F.Rizzi Department of Mechanical Engineering, Johns Hopkins University, Baltimore, MD, USA
R. E. Jones and B. J. Debusschere Sandia National Laboratories, Livermore, CA, USA
O. M. Knio∗ Department of Mechanical Engineering and Materials Science, Duke University, Durham, NC, USA (Dated: April 29, 2013) This article extends the uncertainty quantification (UQ) analysis introduced in Part I for molecular dynamics (MD) simulations of concentration driven ionic flow through a silica nanopore. Attention is now focused on characterizing, for a fixed pore diameter of D = 21 ˚ A, the sensitivity of the system to the Lennard-Jones energy parameters, εN a+ and εCl− , defining the depth of the potential well for the two ions Na+ and Cl− , respectively. A forward propagation analysis is applied to map the uncertainty in these parameters to the MD predictions of the ionic fluxes. Polynomial Chaos (PC) expansions and Bayesian inference are exploited to isolate the effect of the intrinsic noise, stemming from thermal fluctuations of the atoms, and properly quantify the impact of parametric uncertainty on the target MD predictions. A Bayes factor analysis is then used to determine the most suitable regression model to represent the MD noisy data. The study shows that the response surface of the Na+ conductance can be effectively inferred despite the substantial noise level, whereas the noise partially hides the underlying trend in the Cl− conductance data over the studied range. Finally, the dependence of the conductances on the uncertain potential parameters is analyzed in terms of correlations with key bulk transport coefficients, namely viscosity and collective diffusivities, computed using Green-Kubo time correlations.
∗
[email protected]
2 I.
INTRODUCTION
Molecular dynamics (MD) simulations are widely encountered in both industrial and academic environments, playing a key role for the study of a large variety of systems, ranging from liquids to solids, as well as proteins and nucleic acids (DNA, RNA). As every simulation technique, MD is an approximate method, and its key “ingredient” is a potential function that models the interactions between atoms. As a consequence, the reliability of MD critically depends on the accuracy with which the potential function can reproduce the atomic interactions occurring in the real system of interest. Typically, MD potentials are formulated in terms of two-body or many-body contributions, with long-range and short-range terms, each requiring a suitable functional form [1] and specific model parameters. Despite the fact that MD methods have been continuously developed and improved since their introduction in the late 1950s, literature shows that the definition of MD potentials, particularly the choice of the governing parameters, still poses a key challenge, and can represent an important source of uncertainty. A case-in-point is the large variety (more than 50) of MD models proposed in the past trying to suitably describe the water molecule and its interactions in the form of governing potentials, see e.g. the reviews in Refs. [2–5]. A common feature shared by a large subset of these models is that they rely on a Lennard-Jones formulation to describe attraction and repulsion (Van der Waals) forces between the molecules, but the values of its defining parameters vary substantially across the full spectrum of models [2–5]. Most of these potential parameters are obtained by fitting them to reproduce bulk-phase experimental data in classical MD or Monte Carlo simulations. Hence, these potentials are capable of reproducing only some physical properties of water with a good degree of accuracy, while substantial deviations from experiments occur for other properties [3, 5]. Reported experiences also indicate how some of these representations overall have a better performance than others, and, thus, a subset of these models has become widely adopted for computational settings, e.g., the SPC [6] and the TIP4P [7] models. The above discussion grounds the motivation for the present study, which focuses on analyzing how uncertainties in a subset of force-field parameters affect the predicted dynamics for the nanopore MD system introduced in Part I [8] involving a silica pore connecting two reservoirs of an aqueous NaCl solution. A key observation is that, for the present system, the connection between the potential parameters and the observables extracted from the MD simulations is not straightforward, in large part due to the heterogeneous nature of the system and the complexity of the dynamics. This complexity differentiates the present analysis from our prior work [9, 10], which focused on quantifying the effect of potential uncertainties in a pure, homogeneous system. The uncertainty quantification (UQ) analysis presented in this article is based on assuming parametric uncertainty for the Lennard-Jones energy parameters of both ions, thus providing a consistent and self-contained framework. We rely on Polynomial Chaos (PC) expansions [11–13] to represent random variables, and on a Bayesian inference approach [9, 14, 15] to determine the unknown coefficients in the PC expansion. The Bayesian formalism allows us to isolate the effects of thermal noise in the MD computations, and consequently quantify the impact of parametric uncertainty on the MD predictions of the ionic conductances. Attention is then focused on explaining the inferred trends in the ionic conductances in terms of correlations with transport coefficients. The article is organized as follows. In § II, we briefly review the force-field, in § III, we discuss the theoretical framework of the UQ analysis and present the corresponding results in § IV. In § V, we explore the correlation between the observed fluxes and the bulk transport coefficients of the fluid. A brief summary and concluding remarks are finally presented in § VI. II.
REVIEW OF THE POTENTIAL FUNCTION
In this section, we briefly recall the details of the potential formulation, given that it represents the key focus of this second paper. For the details of the MD system and the simulations, the reader is referred to Part I [8]. The total potential energy, Φ, of the system is represented as a combination of pre-defined intra-molecular (bonded ) interactions and inter-molecular (non-bonded ) contributions. For the purposes of this article, we recall only the Van der Waals contribution to non-bonded interactions, [16] which are modeled with a Lennard-Jones (LJ) potential ΦLJ =
NX atoms
4εij
i=1,j>i
|
"
σij rij
12 {z
−
ΦLJ (rij )
σij rij
6 #
}
(1)
where Natoms is the total number of atoms, rij is the distance between the i-th and j-th atoms, ΦLJ (rij ) is the Lennard-Jones (LJ) potential describing attraction and repulsion forces and defined by {εij , σij }, i.e. the depth of the potential well and the distance at which ΦLJ (rij ) becomes zero. As commonly done, we define the LJ parameters
3
Hw Water Ow Ohy Hhy Silica Si Obulk Na+ Ions Cl−
ε [eV] 0.0 0.006721438 0.006595682 0.00199475 0.02602000 0.006595682 0.002033777 0.006504600
σ [˚ A] 0.0 3.15365 3.15380 0.40000 3.91996 3.15380 2.42990 4.04470
Reference Charge |e| Reference [7, 17] 0.520 [7, 17] [7, 17] −1.040 [7, 17] [18] −0.51 [18] [18] 0.32 [18] [18] bond incrementally [19, 20] [18] −0.70 [18] [21] +1 [21] [21] −1 [21]
TABLE I. Lennard-Jones parameters (ε, σ) and Coulombic charge, in multiples of the electron charge |e|, for each atom type present in the system: {Hw , Ow } are water hydrogen and oxygen atoms, respectively, {Hhy , Ohy } are hydrogen and oxygen atoms appearing in an hydroxide group, Obulk are oxygen atoms in the bulk of the silica, Si are silicon atoms and {Na+ , Cl− } are the salt ions.
for the self interactions, and use the Lorentz-Berthelot (LB) mixing rules to calculate the cross-species interaction √ parameters, namely σαβ = 12 (σα + σβ ), and εαβ = εα εβ . For the silica, the corresponding set of potential parameters describing the bonded (internal bonds, angles, dihedrals), and the non-bonded interactions are extracted from Lopes et al. [18]. For water, we assume the TIP4P model [7, 17, 22], and all its potential parameters are extracted from Refs. [7, 17]. Finally, for the ions Na+ and Cl− , which are only subject to non-bonded interactions (LJ and Coulombic), we fix all parameters based on Ref. [21], except for the two LJ energy parameters εN a+ , and εCl− , which will play the role of parametric uncertainty in the present study. Table I summarizes the LJ parameters and the atomic charges for each atom type. For the full details of the force-field the reader is referred to Part I [8]. III.
UQ ANALYSIS
˚, we perform a UQ analysis based on two uncertain potential parameters, namely For a fixed pore diameter D = 21 A εN a+ and εCl− , which define the depths of the LJ potential well for Na+ and Cl− , respectively. This limited choice stems from the large computational cost associated with the present MD system, and is motivated by the anticipated impact of εN a+ and εCl− on the system dynamics, particularly on the MD predictions of the ionic conductance. Our approach is based on a forward propagation of uncertainty, similar to the gating charge analysis discussed in Part I [8]. We parametrize the uncertainty in the target potential parameters, εN a+ and εCl− , by rewriting them as εN a + = µ ε + + α ε + ξ 1 = 0.002033777 + 0.0019923703 ξ1 , εCl− = µε− + αε− ξ2 = 0.006504600 + 0.0008630547 ξ2 ,
[eV],
(2a)
[eV],
(2b)
where {ξ1 , ξ2 } are independent, identically-distributed (i.i.d.) standard random variables uniformly distributed in the interval (−1, 1), the two variables {µε+ , αε+ } identify the mean and standard deviation of εN a+ , while {µε− , αε− } are the mean and standard deviation of εCl− . The reformulation is based on uniform random variables due to the assumption that both εN a+ and εCl− have a compact continuous support. The choice of a uniform density mainly stems from the fact that no a priori information or “expert knowledge” is available about the parameters, and, thus, a uniform density seems to be the most suitable choice, since it is highly non-informative. Alternately, if one had more a priori information about the parameters, for instance, more detailed knowledge of the nature of their uncertainties, one could resort to an more informative distribution, e.g. a Gaussian. The means, µε+ , µε− , are set equal to their reference values used for the study presented in Part I [8], whereas the standard deviations, αε+ , αε− , are calculated based on the values reported in Ref. [21]. The stochastic reformulation in (2) has two main consequences on the force-field. First, it directly affects the depth of the potential well for Na+ –Na+ and Cl− –Cl− LJ interactions. Second, since the cross-interaction parameter, εαβ , between different atom types is calculated using the LB mixing rule, the reformulation in (2) also affects all the cross-interactions between each ion and the other atoms. Figure 1 (a) shows how the uncertainty in (2) affects the form of the Lennard-Jones potential ΦLJ (r), see equation (1), for Na+ (blue) and Cl− (red). Note that the input uncertainty affects only the depth of the potential wells, while the lengths at the which the potentials become zero are fixed to their reference values σN a+ = 2.42990, σCl− = 4.04470. The following observations can be made. First, the uncertain range characterizing εN a+ is substantially larger than that for εCl− , reflecting the current state of knowledge about these MD parameters as shown by
4 x 10
−3
x 10
−3
Range for Na+
2
Range for Cl−
Φ L J (µ ǫ + − α ǫ + )
0 Φ LJ [eV]
Φ LJ [eV]
0 −2
LJ for cross Na+−Cl− interaction
2
Φ L J (µ ǫ + ) Φ L J (µ ǫ − − α ǫ − )
−4 Φ L J (µ ǫ + + α ǫ + )
Φ L J (µ ǫ − )
−6
in m in Φ L J (ǫ m N a +, ǫ C l −)
−2
Nominal curve
−4
ax m ax Φ L J (ǫ m N a +, ǫ C l − )
−6
Φ L J (µ ǫ − + α ǫ − )
−8 1
2
3
4
5
r [˚ A]
6
7
8
−8 1
9
2
3
(a)
4
5
r [˚ A]
6
7
8
9 (b)
FIG. 1. (a) Lennard-Jones potential, ΦLJ (r), computed for Na+ -Na+ (blue) and Cl− -Cl− (red) interactions based on the uncertainty ranges defined for εN a+ , εCl− in equation (2) and for the fixed values σN a+ = 2.42990, σCl− = 4.04470. Panel (b) √ shows the LJ potential, ΦLJ (r), for the cross-interaction Na+ -Cl− where εN a+ Cl− = εN a+ εCl− is calculated using the uncertain ranges in equation (2), and σN a+ Cl− = (σN a+ + σCl− )/2.
Ref. [21]. Secondly, we observe that the LJ-curve obtained for Cl− using the smallest value of εCl− is still deeper than the LJ-curve computed for Na+ using the maximum value of εN a+ . Physically, this means that the strength of the Cl− –Cl− interactions are always larger than those for Na+ –Na+ . As a representative example of a cross-interaction, Figure 1 (b) shows how the uncertainty introduced in (2) shapes the LJ potential governing Na+ -Cl− cross-interaction. This plot is obtained using the LB mixing rule to calculate εN a+ Cl− , and exploiting equation (2) for the uncertainty. Again, the length, σN a+ Cl− , at which the potential becomes zero has no uncertainty and is given by σN a+ Cl− = (σN a+ + σCl− )/2. The plot shows that the mutual interaction between the ions is substantially affected, mainly due to the large uncertainty in εN a+ . Similar plots and analysis can be made for the remaining cross interactions, i.e. between each ion and the remaining atomic types, but, for brevity, are omitted from the present manuscript.
A.
Bayesian regression: formulation
The stochastic reformulation of the inputs in equation (2) implies that the observables extracted from the MD simulations of the nanopore flow can be considered as random variables. If it has finite variance, any quantity of interest (QoI) can be expressed in terms of a suitable polynomial chaos expansion (PCe) [13]. For the present 2D case, the PCe of a generic QoI, X, can be expressed as X=
∞ X
ck Ψk (ξ) ,
(3)
k=0
where {ck }∞ k=0 is the set of PC coefficients, and Ψk are multi-dimensional Legendre basis functions of the “germ” ξ = {ξ1 , ξ2 }, introduced in equation (2). In practice, the basis is truncated at a given polynomial order, p, resulting in an expansion with P + 1 terms, where P + 1 = (2 + p)!/2!/p! [13]. Equation (3) represents a direct relationship between the uncertain inputs εN a+ (ξ1 ) and εCl− (ξ2 ) and a target observable X, and, thus, provides a mapping of the parametric uncertainty from the potential to the observables. To determine the PC coefficients, we recall two main classes of methods. The first class exploits the orthogonality of the Legendre Polynomials to derive an integral expression for each PC coefficient. Examples of this class are the Intrusive Spectral Projection (ISP) and Non-intrusive Spectral Projection (NISP) [12, 13, 23–28]. As in Part I [8], we rely on a Bayesian regression approach [9, 10, 14, 15] to determine the unknown PC coefficients. This involves the following three main steps: collecting a set of the observations {ξ i , Xi } of a target observable, X, by suitably sampling the 2D ξ-space and running the MD system at these sampling points, formulating the Bayesian model and, finally, sampling the target posterior distribution.
1
1
1
0.75
0.75
0.75 0.5 0.25
0 −0.25
ξ2
0.5 0.25 2
0.5 0.25 ξ
ξ
2
5
0 −0.25
0 −0.25
−0.5
−0.5
−0.5
−0.75
−0.75
−0.75
−1 −1 −0.75−0.5−0.25 0 0.25 0.5 0.75 1 ξ
1
−1 −1 −0.75−0.5−0.25 0 0.25 0.5 0.75 1
(a)
ξ1
−1 −1 −0.75−0.5−0.25 0 0.25 0.5 0.75 1
(b)
ξ1
(c)
FIG. 2. (a,b) 2D-Fej´er grids obtained at the second (a) and third (b) resolution levels, showing the nested nature by color-coding red the shared subset of points between levels. Panel (c): sampling grid used in the present study.
1.
Space sampling and collection of observations
To sample the ξ-space, following the approach introduced in Ref. [9], we employ a sampling scheme based on Fej´er grids [13, 29].[30] Figures 2 (a,b) show the 2D Fej´er grids obtained at the second ℓ = 2 (a) and third ℓ = 3 (b) resolution levels. The nested nature makes these grids well suited for successive space refinement as shown in Ref. [9]. In the present study, we sample the ξ-space using the grid shown in Figure 2 (c), which results from the combination of the grid shown in panel (a), with four additional points extracted from that in panel (b), which are used to refine the sampling in the top-left and bottom-right corners. These two opposite corners identify regions of the space where one of the two uncertain LJ parameters (εN a+ , εCl− ) is close to its maximum, while the other is close to its minimum, thus representing regions of particular interest. This grid is a suitable starting point, since an adaptive refinement as shown in Ref. [9], can be further applied, if necessary, by introducing additional nested points from the grid in panel (b). n=13 Each grid point ξ i i=1 corresponds to a particular set of the values of the uncertain LJ parameters as defined by equation (2). To account for the effect of the intrinsic (thermal) noise, we run m = 3 realizations of the MD system at each grid point. These replicas are obtained using 3 different sets of random numbers to initialize the velocity field of the atoms and the insertion/deletion locations and velocities of the concentration control algorithm. Both the sampling grid and the limited number of replicas are chosen to mimic the situation when only a limited number of model runs is available. We could have run more MD simulations to collect more replicas and, thus, reduce a priori the uncertainty, but we prefer to analyze a “small” data set and let the Bayesian framework properly account for the resulting uncertainty. Running the MD simulation for each replica at each sampling point in the ξ-space yields the following set of n × m = 39 observations j=1,...,3 G = Gi,j i=1,...,13 ,
(4)
where, generically, G denotes the steady-state conductance computed for either Na+ or Cl− , the index i enumerates the grid points, j enumerates the replicas, and, consequently, Gi,j denotes the j-th replica value obtained at the i-th sampling point ξ i . Figure 3 shows the data-set, G, obtained for Na+ (a) and Cl− (b), with the observations color-coded based on their location in the ξ-space. The dashed lines are superimposed to visualize the grid point ξ i at which the 3 replica values are obtained. The figures show that the substantial noise level prevents an observer from extracting a distinct trend from the data of both ionic conductances by inspection. A comparison of the plots shows that, for a given sampling node, the magnitude of the Cl− conductance, GCl− , is larger than the corresponding value obtained for GN a+ . This result is expected for the present choice of the pore diameter, D = 21 ˚ A, since in this case the diffusivity of the ions dominates with respect to size (steric) effects, as described in Part I [8]. Therefore, due to the larger diffusivity, the flux of Cl− is larger.
100 90 80 70 60 50 40 30 20 10 0 1
G [nm/ns]
G [nm/ns]
6
0.5
0 ξ2
−0.5
−1 −1
0
−0.5 ξ1
0.5
1 (a)
100 90 80 70 60 50 40 30 20 10 0 1
0.5
0
−0.5
ξ2
−1 −1
0
−0.5 ξ1
0.5
1 (b)
FIG. 3. Data-set, G, obtained for Na+ (a) and Cl− (b), where we color-code the observations based on their location in the ξ-space. The dashed lines are superimposed to visualize the specific grid point ξi , for each of which a set of 3 replica values is obtained. 2.
Likelihood function and Bayes theorem
The regression analysis aims at representing the conductance data, G, with a regression function, M (ξ1 , ξ2 ), which we assume in the form of a Legendre-Uniform PC expansion as M (ξ1 , ξ2 ) =
P X
gk Ψk (ξ1 , ξ2 ).
(5)
k=0
The problem consists in applying Bayesian inference to infer the PC coefficients {gk }P k=0 , given a data-set, G, of the form in equation (4). We assume an additive error regression model according to Gℓ = M (ξ ℓ ) + γℓ ,
ℓ = 1, . . . , 39,
(6)
where the index ℓ enumerates the 39 available data points, ξ ℓ denotes the coordinates of the ℓ-th observation Gℓ , and γℓ is random variable capturing the discrepancy between the data point, Gℓ , and the corresponding model prediction, M (ξ ℓ ). Similarly to the gating charge analysis in Part I [8], since each conductance data point Gℓ is extracted using a running average in independent but statistically equivalent MD simulations, we can leverage central limit arguments to argue that as the number of atoms in the system and the number of time-averaged samples become large, the distribution of the target conductance, G, around the true mean tends to a Gaussian. Based on the previous arguments, we assume 2 {γℓ }N ℓ=1 to be independent and normally distributed with mean zero. To properly model the noise variance, σℓ , of each Gaussian random variable, γℓ , we first need to explore the behavior of the variance extracted from the data. + Figure 4 shows a 3D plot of the variances obtained from the 3 replicas at each grid point {ξ i }13 i=1 for Na (a) and − Cl (b). The plots show that the data seem to outline an approximate linear trend yielding a maximum in the region close to the line ξ1 = 1. To verify whether a linear trend is actually present, we perform a linear fit of the data, and estimate its goodness of fit in terms of the coefficient of determination. This quantity, which ranges between 0 and 1 and has a standard definition widely used in statistics, provides a measure of how well a regression line fits a set of data. A value close to 1 indicates that the regression model fits the data well, while a value close to 0 indicates a negative outcome. For the present case, a linear fit model for the variances of Na+ gives a coefficient of determination of ∼ 0.6, while for the variances of Cl− we obtain ∼ 0.55. When using a cubic model to fit the data, the coefficient of determination for the Na+ variances rises to ∼ 0.7 and that for Cl− to ∼ 0.75, but at the expense of having 7 more coefficients to include in the analysis. These results show that a linear trend represent a good compromise between accuracy and simplicity, and, thus, suggest to use in the inference a linear ξ-dependent noise model. To this end, we parametrize the standard deviation, σ(ξ), with a linear PC expansion of the form σ(ξ) = h0 + h1 Ψ1 (ξ) + h2 Ψ2 (ξ),
(7)
where {hk }2k=0 are the PC coefficients, which are treated as hyperparameters. As will be shown, and in light of the preceding regression analysis performed, this parametrization is sufficient to provide a good representation of the observed trends in the variance.
150 135 120 105 90 75 60 45 30 15 0 1
Variance [(nm/ns) 2 ]
Variance [(nm/ns) 2 ]
7
0.5
0 −0.5 −1 −1 ξ2
−0.5 0 ξ1
0.5
1 (a)
150 135 120 105 90 75 60 45 30 15 0 1
0.5
0 −0.5 −1 −1 ξ2
−0.5 0 ξ1
0.5
1 (b)
+ FIG. 4. 3D plot of the variances obtained from the 3 replicas run at each grid point {ξi }13 i=1 . The results are shown for Na (a) and Cl− (b), and are color-coded based on their location in the ξ-space. The dashed lines are superimposed to visualize the grid point ξi at which the variance is computed.
Again, as done in Part I [8], we omit a model discrepancy term in the error formulation (6), because we prefer to keep it as simple as possible, and since we consider the noise variance as a hyperparameter. This assumption reflects the fact that, for a given observable, the regression function, M (ξ), is assumed to properly capture the corresponding data set, G. This hypothesis can be verified a posteriori by comparing the trend and magnitude inferred for σ 2 with its estimate based on the data. For a given order, p, of the PC regression function, M (ξ), if the inferred behavior of σ 2 (ξ) is in good agreement with its corresponding data-based estimate shown in Figure 4, it follows that the model representation is appropriate. On the contrary, if significant differences arise, the regression function must be refined in complexity since a model discrepancy effect between the approximation and the true behavior is lumped into the data noise. This line of reasoning leads to the following likelihood function 2 P gk k=0 , hl l=0 ; G 3 13 Y Y 1 [Gi,j − M (ξ i )]2 p = , exp − 2[σ(ξ i )]2 2π[σ(ξ i )]2 i=1 j=1
L=L
(8)
and Bayes’ theorem then yields the joint posterior
P 2 Y 2 2 Y P P π gk k=0 , hl l=0 G ∝ L gk k=0 , hl l=0 ; G qk g k qˆl hl , k=0
(9)
l=0
where qk gk and qˆl hl denote, respectively, the presumed independent priors of the regression function PC coeffi P 2 cients, gk k=0 , and of the coefficients, hl l=0 , defining the PCe of the noise standard deviation. For all priors, we assume uniform distributions with suitably large bounds. Moreover, since we adopted a linear PCe for the noise standard deviation, σ(ξ), but the variance σ 2 (ξ) is the quantity appearing in the likelihood, we properly constrain the priors, {ˆ ql }2l=0 , of {hl }2l=0 to be either all positive or all negative, thus avoiding the possibility of obtaining a multimodal distribution. We remark, however, that since in this case no prior information about the coefficients is available, this choice is arbitrary. For verification purposes, we also ran the Bayesian inference using Jeffreys priors[31] for all target parameters, and this approach led to similar results. For brevity, we omit this discussion from the present manuscript. The posterior density in equation (9) is sampled with a Markov chain Monte Carlo (MCMC) method based on the adaptive Metropolis (AM) algorithm [32–36], using the same setting described in Part I [8]. Each MCMC chain is run for 35000 steps, and the first 15000 steps are discarded to eliminate the “burn-in” period. For brevity, we omit the plots and discussion of the MCMC chains obtained during the inference, but confirm that in all cases, the runs are repeated and tuned to ensure that all the chains are well-mixed.
8 IV.
UQ RESULTS
For each conductance, the inference analysis is run with a linear PCe for the noise standard deviation, as shown in equation (7), but exploring four subsequent orders, p, of the regression function, M (ξ). We consider a constant (p = 0), linear (p = 1), quadratic (p = 2), and cubic (p = 3) PC model. Note that higher-order (p ≥ 4) PC models cannot be employed since the problem would be under-determined, because the number (n = 13) of sampling nodes in the ξ-space would not be sufficient to cover the corresponding degrees of freedom. 0.08
0.08
p=0 p=1 p=2 p=3
0.07 0.06
0.06 0.05 π (h 20 )
π (h 20 )
0.05 0.04
Data−based estimate
0.03
0.02
0.02
0.01
0.01 25 50 75 100 125 150 175 200 225 250 275 300 (h 0 ) 2 [(nm/ns) 2 ] (a)
Data−based estimate
0.04
0.03
0 0
p=0 p=1 p=2 p=3
0.07
0 0
25 50 75 100 125 150 175 200 225 250 275 300 (h 0 ) 2 [(nm/ns) 2 ] (b)
FIG. 5. Marginalized posterior, π(h20 ), obtained via KDE from the chain samples of the squared mean (h0 )2 of the PCe of the noise standard deviation, computed for Na+ (a) and Cl− (b). The results are shown for inferences based on using four different orders p = 0, 1, 2, 3 (black,red,green,blue) of the target regression function, M (ξ).
As a first representative result, Figures 5 (a,b) show the marginalized posterior, π(h20 ), obtained via Kernel Density Estimation (KDE) using the chain samples obtained for the squared mean of the PCe of the noise standard deviation. The results are shown for Na+ (a) and Cl− (b), and are based on four different PC orders p = 0, 1, 2, 3. A dashed line in the figure indicates the data-based estimate of the variance which, in this case, is set equal to the mean of the variances obtained at the 13 grid points in the stochastic space. Figures 5 (a,b) allow us to compare the inferred value of the noise variance, σ 2 , with its estimate based on the data and, thus, draw some preliminary conclusions regarding the accuracy of the regression function considered in the inference. The results computed for Na+ , see Figure 5 (a), show that the marginalized posteriors, π(h20 ), obtained for p ≥ 1 are peaked nearly exactly around the data-based estimate of the variance. On the contrary, the posterior, π(h20 ), obtained using a zeroth-order (p = 0) regression function is substantially shifted toward larger values. This indicates that a constant-value regression function, M (ξ), does not yield an accurate representation of the data, and so the corresponding discrepancy between M (ξ) and the conductance data is lumped into the data noise, raising its inferred magnitude. These results suggest that an underlying trend is present in the data set obtained for GN a+ , which can be captured with sufficient accuracy using a regression function of order p ≥ 1. For the Cl− conductance, GCl− , Figure 5 (b) shows that regardless of the order, p, of M (ξ), the mode of the posterior π(h20 ) captures the data-based estimate of the variance with good accuracy. Hence, this result indicates that the data-set obtained for GCl− does not lead to a net underlying trend, since a constant value seems to be a suitable representation. The main advantage of a Bayesian inference approach is that it leads to PC coefficients that are random variables characterized by a well-defined joint probability distribution. It follows that this posterior can be exploited to extract important information about its shape and, also, potential correlations between the PC coefficients. To this end, we use the MCMC samples of the third-order (p = 3) PC spectrum obtained for Na+ and Cl− to compute the matrix of correlation coefficients, ζ(gi , gk ), i, k = 0, . . . , 9 for the PC spectra. The correlation coefficient, ζ(X, Y ), of two random variables X and Y is defined as ζ(X, Y ) =
E[(X − µX )(Y − µY )] , σX σY
(10)
where E is the expectation operator, while (µX , σX ) and (µY , σY ) are means and standard deviations of X and Y , respectively. It falls in the range (−1, 1), with 1 indicating a perfect correlation and −1 for a case of perfect anticorrelation, and its sign determines the tendency in the linear relationship between the variables. The top row of
9 Figure 6 shows the results of ζ(gi , gk ), i, k = 0, . . . , 9, obtained for Na+ (a) and Cl− (b), with each matrix entry colorcoded based on the corresponding value of ζ(gi , gj ). The following observations can be made. First, the correlation structure arising in the results for Na+ is very similar to the one obtained for Cl− . This is particularly evident for the PC coefficients with indices k > 5, while a more substantial difference is observed for the low-order coefficients (i.e. k < 5). This indicates that the low order coefficients play a key role in differentiating the trends underlying the Na+ and Cl− data. Secondly, for both ions the results show that most correlation coefficients are positive, thus indicating that an increasing trend of gi is (generally) associated with an increasing trend in gj . Negative correlation coefficients are mainly localized in the upper right corner of the plots, see Figures 6 (a,b), which involve the highest-order PC coefficients. Finally, we observe that the strength of the correlations are overall relatively weak, since the magnitude of ζ(gi , gj ) is mainly contained in (−0.5, 0.5) with only few values exceeding this range. The bottom panels in Figure 6 show few sample contour plots of the joint posterior densities π(gi , gk ), i = 0, . . . , 4, k = 1, . . . , 5, using the results inferred for Na+ . Note that these plots correspond to the combinations of coefficients highlighted in Figure 6 (a) with a black dashed lined. The correlation structure obtained for Cl− (not shown) was found to be similar to that of Na+ . For visualization purposes, the range adopted for the y-axis is the same in all plots, while the one for the x-axis varies based on the order but maintains the same width in all plots. This setting allows us to have a clear representation of the shape of these posteriors, and, at the same time, to visually identify the differences in the correlations, and readily make the reference with Figure 6 (a). The joint posteriors reveal, in all cases, profiles closely resembling Gaussians, and evident differences in the underlying correlation structures. In all cases, consistent with the preceding discussion, the joint posteriors display weak correlations. These considerations allow us to conclude that correlations are present in the inferred PC spectra for both ions, but they are overall weak, and the shapes of the joint posteriors, in most cases, can be suitably approximated with Gaussians. It follows that, for the present case, one could resort to using the marginalized posteriors for extracting statistics rather than keeping the full P + 1-dimensional joint density. A.
Bayes factor and model selection
To select the order of the regression function most suitable to describe the MD observations we resort to a model selection analysis using Bayes factors [37]. The Bayes factor, already employed in Part I [8], is a non-dimensional quantity allowing one to compare two “models” describing a given set of data, providing a quantitative measure to discriminate between the two. We recall that the term “model” refers to the complete set of parameters that are being inferred. In the present case, a “model” refers to the parameter vector θ p = {g0 , . . . , gP , h0 , h1 , h2 }, which comprises the set of P coefficients {g0 , . . . , gP } defining the regression function, M (ξ), for a given expansion order p, and the three PC coefficients {h0 , h1 , h2 } defining the linear PCe of the noise standard deviation. Given two different models, namely θ p1 and θ p2 , associated with regression functions of order p1 and p2 , respectively, the Bayes factor, B(θ p1 , θ p2 ), is given by R L (θ p1 ; G) Pr(θ p1 ) dθ p1 B(θ p1 , θ p2 ) = R , (11) L (θ p2 ; G) Pr(θ p2 ) dθ p2
where G is the conductance data-set used in the inference, L (θ ; G) and Pr(θ) are the likelihood function and prior, respectively, previously discussed in equation (9). The Bayes factor has to be computed numerically and, for convenience, we report and discuss its natural logarithm. The more positive the value of ln(B(θ p1 , θ p2 )), the stronger the evidence supporting θ p1 , while the converse applies for negative values [37]. The values of ln(B(θ p1 , θ p2 )), computed for a given ion and all combinations of orders in the range 0 ≤ p1 , p2 ≤ 3 are shown in Table II. Due to the use of the natural logarithm, the results obtained for a given ion yield an antisymmetric matrix of values. The results obtained for Na+ indicate that a second-order (p = 2) regression function yields the most suitable regression model. However, while the evidence supporting this model is decisive when compared to those based on a constant (p = 0) and a cubic (p = 3) expansion of M , it is only minimally dominating with respect
TABLE II. Computed values of ln B(θ p1 , θ p2 ) for all possible combinations of the orders {p1 , p2 } = 0, 1, 2, 3. p1 p1 p1 p1
p2 = 0 =0 – = 1 19.341 = 2 19.933 = 3 16.285
N a+ p2 = 1 p2 = 2 -19.341 -19.933 – -0.593 0.593 – -3.055 -3.648
p2 = 3 -16.285 3.055 3.648 –
p2 = 0 – -2.284 -2.737 -6.499
Cl− p2 = 1 p2 = 2 2.484 2.737 – 0.254 -0.254 – -4.015 -3.761
p2 = 3 6.499 4.015 3.761 –
10 1
g9 g8
g8
g7
g7
0.5
g6
|
|
g5
0
g4
|
|
g3
g2
g2
−0.5
g1 g0
g1
g2
Linear
|
g3
g4
g5
Quadratic
|
g0
g6
g7
g8
g0
−1
g9
g0
Cubic
g1
(a)
g1
g2
Linear
g2
|
g3
g4
g5
Quadratic
|
g6
g3
15
15
15
15
10
10
10
10
5
5
5
5
5
0
0
0
0
0
−5
−5
−5
−5
−5
−10
−10
−10
−10
−10
−15
−15
−15
−15
40
45
50
0
5
10
15
20
−10
−5
0
5
10
−10
15
15
15
15
10
10
10
10
5
5
5
5
0
0
0
0
−5
−5
−5
−5
−10
−10
−10
−10
−15
−15
−15
30
35
40
45
50
0
5
10
15
20
−10
15
15
15
10
10
10
5
5
5
0
0
0
−5
−5
−5
−10
−10
−10
−15
−15
30
35
40
45
50
0
15
15
10
10
5
5
0
0
−5
−5
−10
−1
g9
(b)
−15 −5
0
5
10
−5
0
5
10
−10
−5
0
5
10
−15 −5
0
5
10
−5
0
5
10
−10
−15 5
10
15
20
5
10
15
20
−10
−10
−15 30
g8
g4
10
35
g7
Cubic
15
30
g2
−0.5
g1
g0
g3
0
g4
g3
g4
0.5
g6
g5
g5
1
g9
−15 35
40
45
50
35
40
45
50
0
15 10 5
g1
0 −5 −10 −15 30
FIG. 6. Top row: matrix of correlation coefficients, −1 ≤ ζ(gi , gj ) ≤ 1, i, j = 0, . . . , 9, computed for the third-order PC spectra inferred for Na+ (a) and Cl− (b). Each matrix entry is color-coded based on the corresponding value ζ(i, j). Since the matrix ζ(gi , gj ) is symmetric, only the upper triangular part is shown. Bottom rows: sample joint posterior densities π(gi , gj ), i = 0, . . . , 4, j = 1, . . . , 5, computed for the PC coefficients inferred for Na+ , corresponding to the cases highlighted in panel (a) by a black dashed lined.
the linear case. This indicates that even a linear regression function, M , would be a suitable choice for this case. For the negative ion, Cl− , the results suggest a constant (p = 0) regression function as the most suitable model, which strongly dominates with respect to all other orders. This analysis allows us to conclude that a quadratic (p = 2) regression function would be appropriate to represent the data set of the Na+ conductance, GN a+ , while a constant value (p = 0) should be adopted to describe the data-set
11 for GCl− . We remark that these conclusions result from integrating over the complete parameter space involved in the inference. Hence, they not only reflect the accuracy with which the PC regression function, M (ξ), describes the data points of the conductance, G, but also the accuracy of the inferred representation of the noise standard deviation.
B.
Model uncertainty
The Bayesian formalism yields an uncertain PC representation, M , since its PC coefficients are random variables defined by a joint posterior distribution π({gk }P k=0 ). We now focus on analyzing, for both conductances, the posterior uncertainty in the corresponding inferred regression function M . As discussed, a quadratic (p = 2) representation for M should be adopted to describe GN a+ , while a constant model (p = 0) should be used to represent GCl− . For the sake of the analysis, however, we employ a quadratic expansion for both GN a+ and GCl− , which will provide additional support to the results of the Bayes factor. The top row of Figure 7 shows the marginalized posteriors, π(gk ), obtained via KDE for a quadratic expansion, M , for Na+ (a) and Cl− (b). The bottom row of Figure 7 reports key statistics for the corresponding PC spectra, {gk }P k=0 , showing the means (black circle) and the limits ±2˜ σ (error bars), where the mean and the standard deviation, σ ˜ , are both extracted from the posteriors in Figures 7 (a,b). Panels (a,b) show that for both ions, the posterior distribution of each PC coefficient is fairly symmetric, closely resembling a Gaussian, further confirming the results previously obtained from the joint posterior densities. The figures also show that as the order of the coefficient increases, i.e. from the leading, g0 , to the last coefficient, g5 , the peak of the corresponding posterior shifts towards smaller values, finally settling around zero. This illustrates the rapid decay of the PC spectrum, which is clearly illustrated in Figures 7 (c,d). The plots show that the decay of the PC coefficients magnitude is relatively gradual for the Na+ results, while being substantially sharper for the Cl− results, see Figure 7 (d). More specifically, for the Na+ results, Figure 7 (c) shows that the first two coefficients, g0 and g1 , dominate the spectrum, since the higher-order coefficients oscillate around zero. This observation explains why the previous Bayes factor analysis showed a minimal difference between a linear and a quadratic expansion for the GN a+ , see left column of Table II. A similar consideration can be made for the Cl− results shown in Figure 7 (d). In this case, the spectral decay is even more evident, since g0 is strongly dominating in magnitude, while all the remaining coefficients gk ≥ g1 are considerably smaller in magnitude. This is consistent with the Bayes factor analysis, which showed that a constant value is the proper order to use for the GCl− , see right column of Table II. Finally, note that Figures 7 (c,d) show, for both ions, that the uncertainty in the coefficients increases for the higherorder coefficients, while being relatively small for the linear terms. A comparison of Figures 7 (c,d) demonstrates that the posterior uncertainty obtained for the Cl− coefficients is larger than for Na+ . This is in accordance with the variances plotted in Figure 4 showing that the GCl− conductance data have larger spread than those for GN a+ , and, also, in agreement with the results illustrated in Part I [8], where we showed that Cl− fluxes yield larger noise than that for Na+ . Similar observations hold for the third-order PC coefficients; discussion of the corresponding results is consequently omitted.
1.
Model prediction uncertainty
Figure 8 (a) shows 500 second-order (p = 2) sample spectra, gi = {g0,i , . . . , g5,i }500 i=1 , obtained for the PC regression function, M , inferred for GN a+ , by sampling the corresponding marginalized joint posterior π(g). Superimposed to the plots are the means (black circles) and error bars for the limits ±2˜ σ , extracted from Figures 4 (c,d). Figure 8 (a) reflects the observation that the posterior uncertainty in the PC spectrum is larger for the higher-order terms. This is reflected in the plot by the lines of the spectra spreading out in correspondence to the larger modes. As observed earlier, however, the spectrum shows a sharp decay, and thus, we expect the response surface, M , to be mainly dominated by the relatively small uncertainty in the linear coefficients. To assess the accuracy with which the inferred regression function, M , can represent the original noisy data, we extract 5 sample spectra from those shown in Figure 8 (a), and use them to compute the corresponding model predictions over the stochastic space. The results are reported in Figure 8 (b), where different colors are used to distinguish the 5 response surfaces thus obtained, and the noisy data GN a+ (black circles) are superimposed for reference. The plot shows that all 5 response surfaces capture well the noisy data, but reveal evident differences between one another, due to the substantial noise present in the data. A key observation in Figure 8 (b) concerns the evidence of an underlying trend that GN a+ increases as εN a+ (i.e. ξ1 ) increases. This trend only emerges at this point because the noise present in the original data (Figure 3) partially hides its presence. Figure 8 (b) also reveals that the response surface is weakly sensitive to variations along ξ2 , suggesting that GN a+ is weakly sensitive to changes in εCl− over the range studied.
k=3
40
gk [nm/ns]
gk [nm/ns]
80 75 70 65 60 55 50 45 40 35 30 25 20 15 10 5 0 −5 −10 −15 −20
Linear
0
1
k=4
k=5
π(gk )
0.6 k=0 k=1 k=2 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 −20 −10 0 10 20 30
50
60
70
80
Quadratic
2
3
Index, k
0.6 k=0 k=1 k=2 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 −20 −10 0 10 20 30
4
5 (c)
k=3
40
gk [nm/ns]
(a)
gk [nm/ns]
π(gk )
12
80 75 70 65 60 55 50 45 40 35 30 25 20 15 10 5 0 −5 −10 −15 −20
Linear
0
1
k=4
50
60
k=5
70
80 (b)
Quadratic
2
3
Index, k
4
5 (d)
FIG. 7. Top row: marginalized posteriors, π(gk ), of the PC coefficients, gk , of the regression function M obtained for a quadratic expansion (p = 2) for both Na+ (a) and Cl− . Bottom row: errors bars for the PC spectra, {gk }P k=0 , showing the means (black circle) and the limits ±2˜ σ , where both the mean and the σ ˜ are computed from the posteriors in panels (a,b).
To further investigate the results, Figure 9 shows the contour (a) and 3D (b) plot of the quadratic response surface, M , inferred for GN a+ , evaluated using the maximum a posteriori (MAP) estimates of the PC coefficients. The MAP estimates are chosen as representative values, but we verified that a similar results is obtained using the posterior means. This stems from the fact that the posteriors densities of the PC coefficients are relatively symmetric, as shown before, implying that means and MAP estimates nearly coincide. Figure 9 highlights the trend found for the response surface of the Na+ conductance, GN a+ , showing that it reaches its maximum near the bottom-right corner of the space, i.e. where εN a+ reaches its maximum, while εCl− is at its minimum. In the bottom-left corner, i.e. where both εN a+ and εCl− reach their minima, the predicted values of GN a+ are small. The contour lines also suggest that in the center of the domain, i.e. in correspondence to the “nominal” values of εN a+ , the response surface suggests that GN a+ is weakly sensitive to εCl− . This sensitivity amplifies near the corners, where small variations of ξ2 yield sharper changes in the model response. The preceding analysis can be repeated for the Cl− -conductance, GCl− . In this case, however, the results are not as significant, since a constant regression model is found to be most suitable. For brevity, we only report in Figure 10 the model predictions obtained by evaluating the zeroth-order PCe M (ξ) for 5 samples of the PC coefficient drawn from the posterior π(g0 ). The noisy data (black circles) of GCl− are superimposed for reference. The constant-value response surfaces obtained for GCl− capture well the noisy data, and display, overall, small differences, reflecting the small posterior uncertainty affecting g0 , previously shown in Figures 7 (d).
60 55 50 45 40 35 30 25 20 15 10 5 0 −5 −10 −15 −20
100 M (ξ 1 , ξ 2 ) and GN a + [nm/ns]
g k [nm/ns]
13
80 60 40 20 0
0
1
2
3
4
Index, k
1
5
0.5
0 −0.5 −1
ξ
−0.5
2
(a)
0
ξ
0.5
1
1
(b)
FIG. 8. Panel (a) shows 500 sample spectra, gi = {g0,i , . . . , g5,i }500 i=1 , of the regression function M inferred for GN a+ , obtained by sampling the corresponding marginalized joint posterior π(g). The error bars denote the mean and ±2˜ σ . Panel (b) shows the 3D plot of 5 sample response surfaces M obtained by evaluating M for 5 sample-spectra extracted from those in panel (a). The noisy data used in the inference, GN a+ , are plotted with black markers for comparison. Units [nm/ns]
1
M AP (M ) and GN a + [nm/ns]
55 50
0.5
ξ
2
45 0 40 −0.5
35
55
100 80
50
60
45
40
40
20
35
0 −1 −1
30 −0.5
0
ξ1
0.5
1 (a)
1 0.5 0 −0.5 −1 −0.5 ξ 2
0
ξ
1
0.5
1
30 (b)
FIG. 9. Contour (a) and 3D plot (b) of the response surface, M , inferred for GN a+ based on the MAP estimates of the PC coefficients. In panel (a), the black circles indicate for reference the location of the sampling nodes; in panel (b) the data points used for the inference are plotted with black circles.
V.
BULK TRANSPORT COEFFICIENTS AND CORRELATIONS
˚) pore-fluid system, i.e. a pore large enough such that the flow is For the present non-sterically limited (D = 21 A not limited by the size of the ions as described in Part I [8], a direct connection between the physical effects of LJ parameters, εN a+ and εCl− , on the pore flow are difficult to elucidate directly. This is the case because whereas the role of εN a+ and εCl− is clear at the particle interactions level, the conductance, G, is a large-scale observable which results from the non-trivial combination of several intermediate effects. Consequently, in this section we resort to exploring the impact of εN a+ and εCl− on bulk transport properties, namely viscosity and diffusion. Our goals are twofold: first, to investigate the sensitivity of the bulk transport coefficients to the uncertain force-field parameters; and second, to explore, for both ions, the potential correlation existing between the the estimated transport coefficients and the conductances. The bulk transport coefficients for the NaCl solution are computed in a separate MD study. We consider a cubic domain of dimensions 403 ˚ A3 , with periodic boundary conditions along each axis. In order to study the same fluid
14
M (ξ 1 , ξ 2 ) and GC l − [nm/ns]
100 80 60 40 20 0 1
0.5
0 −0.5 −1
ξ2
−0.5
0
ξ1
0.5
1
FIG. 10. 3D plot of 5 zeroth-order sample response surfaces of M obtained for the Cl− conductance by using 5 samples extracted from the posterior π(g0 ). The superimposed noisy data used in the inference, GCl− , are depicted with black markers for comparison.
used to run the pore flow analysis, the cubic domain is filled with an aqueous, 1.5 mol/l NaCl solution, and the TIP4P model is adopted for water. The system comprises a total of 6351 atoms. We focus on analyzing the viscosity of the solution, and the collective diffusion coefficients.[38] In the dilute limit, both transport properties are directly related to the ion mobility: the viscosity is inversely proportional, while the diffusion is directly proportional to the mobility. The target transport coefficients are computed using the Green-Kubo equations, outlined in Appendix A. For notational convenience, we denote µ as the dynamic viscosity, whereas the collective diffusivities are expressed as D++ , D−− and D+− . The first two, D++ and D−− , represent, respectively, the diffusivity of the positive, Na+ , and negative, Cl− , ions due to concentration gradients for the same species. The last quantity, D+− , represents the diffusivity of Na+ due to cross-interactions with Cl− , and vice versa. The symmetry of the collective diffusion matrix, i.e. D+− = D−+ , is an Onsager reciprocal relationship and is reflected in the formulation adopted here to compute the coefficients, see Appendix A. For the force-field, we adopt the same setting governing the NaCl-solution in the pore flow simulation in § II. To ensure the same UQ setting adopted for the analysis of the conductance, we explore the sensitivity of the bulk properties to εN a+ and εCl− , with the same stochastic reformulation introduced in equation (2). The goal is to use Bayesian inference to obtain a PC representation for the bulk transport coefficients as a function of the two uncertain potential parameters. A.
Bayesian regression
The analysis is similar to the one previously described for the pore flow system. The stochastic space, (ξ1 , ξ2 ), is sampled with the same grid employed for the analysis of the conductance shown in Figure 2 (c). At each grid point in the stochastic space, we run 10 independent replica simulations of 1 ns each, using a time step of 0.001 ps, and, using the techniques described in the Appendix, we compute the four transport coefficients of interest. The choice of using 10 replicas is due to the fact that these MD runs are not as expensive as the nanopore flow, and allows us to better account for noise. We thus obtain for each transport coefficient, a set of 130 observations, which can be compactly written as
++ −− +− j=1,...,10 , µ, D++ , D−− , D+− = µi,j , Di,j , Di,j , Di,j i=1,...,13
(12)
where the index i enumerates the n = 13 grid points in the stochastic domain, see Figure 2 (c), j enumerates the 10 replicas, and, consequently, the subscript i,j denotes the j-th replica value of a given observable obtained at the i-th sampling point ξ i = (ξ1,i , ξ2,i ). Each data set is exploited with a Bayesian regression analysis to obtain an individual PC expansion for each transport coefficient. For generality, let F (ξ1 , ξ2 ) be the PC model of a transport coefficient of interest. The PC representation of F is given by PF
. X fk Ψk (ξ1 , ξ2 ), F (ξ1 , ξ2 ) = k=0
(13)
15 where fk if the k-th PC coefficient, Ψk is the k-th basis function and PF denotes the limit of the expansion for a target order, pF , of the basis functions. The four sets of PC coefficients for the transport coefficients are inferred with a Bayesian regression approach, similar to the one applied to the pore conductance in § III A. The only difference is that we now adopt a simplified approach that involves a noise model independent of ξ-space, which can be viewed as a special case of the model in § III A corresponding to using h1 = h2 = 0 in equation (7). This choice is motivated by the fact that the variance computed for the replica values of the transport coefficients undergoes minimal variations over the stochastic space, and, thus, a constant value assumption is a suitable choice. To avoid redundancy, we omit the theoretical details from the present discussion, and directly illustrate the results.
B.
Model selection and results
For each transport coefficient, in order to select which “model” is most suitable to represent the data and the associated information, we again rely on model selection via Bayes factor [37]. For brevity, we limit the discussion to the results. The analysis yields that a linear (p = 1) expansion appears to be the best choice for the viscosity, µ; a quadratic (p = 2) PC model is suitable to represent both the Na+ diffusivity, D++ , and the cross-diffusion coefficient, D+− ; finally, a third-order (p = 3) expansion is the most appropriate for the Cl− diffusivity, D−− . Figure 11 shows the 3D plots of the response surface, F , inferred for µ (a), D++ (b), D−− (c) and D+− (d), computed using the MAP estimates of the corresponding inferred PC coefficients. Figure 11 (a) shows that µ reaches its maximum near the point (ξ1 , ξ2 ) = (−1, −1), indicating that it increases as the LJ parameters εN a+ and εCl− decrease. It reaches its minimum at (ξ1 , ξ2 ) = (1, 1), indicating that it is small for the maxima of εN a+ and εCl− . Figure 11 (b) shows that the Na+ diffusivity is weakly sensitive to variations in ξ2 , and is mainly affected by changes along ξ1 . Specifically, the figure shows that D++ reaches its largest value near the region ξ1 = 1, and its minimum for ξ1 = −1. This indicates that D++ is mostly sensitive to variations in εN a+ and minimally affect by εCl− over the range studied. Figure 11 (c) shows that the response surface D−− , that is described by a cubic PCe, is characterized by being fairly constant throughout the space, with a sharp drop near the (−1, −1) corner. Finally, Figure 11 (d) shows that the mutual diffusivity is mainly sensitive to variation in ξ2 , and is large in the central region of the space, while decreasing near the the regions for ξ2 = 1 and ξ2 = −1. C.
Analysis of conductance-transport correlation
We now focus on exploring the correlations arising between the conductances and transport coefficients to explain the trends previously observed in the conductances in terms of the behavior of the transport coefficients over the stochastic space. The focus of this section is to correlate each conductance with each of the four transport coefficients. This can be achieved in an elegant way by exploiting the PC representation of each observable previously inferred with Bayesian regression. Recall that we used M (ξ) (see equation 5) to denote the PC model inferred to represent the Na+ or Cl− conductance, and F (ξ) (see equation 13) to denote the PC representation of a target transport coefficient. These PC expansions can be exploited to express the covariance, Cov(M, F ),[39] between an individual conductance, M , and an individual transport coefficient, F , as Cov(M, F ) = E (M − E[M ])(F − E[F ]) " P X # PF X =E gk Ψk (ξ) fl Ψl (ξ) k=1
=
PF P X X
k=1 l=1
E gk fl Ψk (ξ)Ψl (ξ)
min(P,PF )
=
X
k=1
l=1
gk fk E Ψk (ξ)2 ,
(14)
where the first step results from the fact that, by definition, the expectation of a PC expansion equals the zeroth-order term, i.e. E[M ] = g0 and E[F ] = f0 , the second step stems from the linearity of the expectation operator, and, the final step stems from the orthogonality of the basis functions. Equation (14) shows that the covariance can be interpreted as a linear combination of the product of the PC coefficients, where the norms of the PC basis functions, E Ψk (ξ)2 , play the role of the weights.
16 −3
x 10
−4
−3
1.5
9.8 9.75
1.25
9.7
1
9.65 0.75
9.6
0.5 1 0.5
9.55 9.5 0 −0.5
ξ2
−1 −1
−0.5 0
0.5
1
ξ1
M AP (F ) and D−− [mm 2 /s]
5 4
0.004 0.002
3
0 1 0.5 0 −0.5 ξ −1 −1 2
2 −0.5 0 ξ1
0.006
2.5
0.004 0.002 2
0 1 0.5 0 −0.5 ξ −1 −1 2
−0.5 0
0.5 1 1.5
ξ1
(b)
x 10
0.01 0.006
0.008
−3
6
0.008
3
0.01
(a)
x 10 0.012
0.012
9.45
M AP (F ) and D+− [mm 2 /s]
M AP (F ) and µ [Pa · s]
x 10
M AP (F ) and D++ [mm 2 /s]
x 10 9.85
0.5 1 1 (c)
−3
7
0.012 0.01
6.5
0.008 0.006
6
0.004 5.5
0.002 0 1 0.5 0 −0.5 ξ −1 −1 2
5
−0.5 0 ξ1
0.5 1 4.5 (d)
FIG. 11. 3D plots of the response surface, F , obtained for µ (a), D++ (b), D−− (c) and D+− (d), computed using the MAP estimates of the corresponding PC coefficients. Each plot shows the 10 × 13 data points (black circles) used for the inference, while the dashed lines link each set of 10 replicas to the corresponding sampling node in the stochastic space.
In the present work, the PC representations are uncertain, hence the covariance in equation (14) between any pair of observables is a random variable. To properly characterize this uncertainty, we adopt the following sampling approach: first, given a pair of observables, M and F , we can draw two sample PC spectra, gi and fi , from the corresponding posteriors π(g) and π(f ); second, each sample spectrum, gi and fi , can be used to evaluate the corresponding covariance with equation (14). This process can be repeated for as many times as desired. The resulting uncertainty in the covariance thus reflects the posterior uncertainty in the PC models which, in turns, reflects the effect of the intrinsic noise. Given that the covariance is a dimensional quantity, it is preferable to consider the correlation coefficient Θ(M, F ) = Cov(M, F )/σM σF , where σM and σF are the standard deviations calculated for the PCe’s M and F , respectively. For the sake of notation, we denote with M − the constant-value PC representation inferred for the Cl− conductance, GCl− ; with M + we denote the quadratic PC representation inferred for GN a+ ; with F µ the linear PCe inferred for µ; with F ++ the quadratic PCe inferred for D++ ; with F −− the cubic PCe inferred for D−− and, finally, with F +− the quadratic model inferred for D+− .
5000
5000
4000
4000
3000
3000
frequency
frequency
17
2000 1000
1000
0 −1 −0.8 −0.6 −0.4 −0.2
0
0.2 0.4 0.6 0.8
Θ(M + , F µ )
0 −1 −0.8 −0.6 −0.4 −0.2
1
5000
5000
4000
4000
3000
3000
2000 1000
0
0.2 0.4 0.6 0.8
Θ(M + , F ++)
(a)
frequency
frequency
2000
1 (b)
2000 1000
0 −1 −0.8 −0.6 −0.4 −0.2
0
+
0.2 0.4 0.6 0.8
Θ(M , F
)
−−
0 −1 −0.8 −0.6 −0.4 −0.2
1 (c)
0
0.2 0.4 0.6 0.8
Θ(M + , F +−)
1 (d)
FIG. 12. Histograms based on the 50000 samples obtained for Θ(M + , F µ ) (a), Θ(M + , F ++ ) (b), Θ(M + , F −− ) (c), and Θ(M + , F +− ) (d).
D.
Correlation between M − and transport coefficients
Due to the fact that the Cl− -conductance, GCl− , is represented with a constant-value PC representation, M − , its correlation coefficient computed with each transport coefficient is zero, i.e. Θ(M − , F µ ) = Θ(M − , F ++ ) = Θ(M − , F −− ) = Θ(M − , F +− ) = 0. This inconclusive result stems from the fact that a trivial representation of the Cl− conductance has been obtained as a result of the effect of the intrinsic noise, as discussed in § III A. E.
Correlation between M + and the transport coefficients
To provide a suitable characterization of the uncertainty in the covariance, we proceed as follows. First, using the joint posterior distribution, π(g0 , . . . , g5 ), of the quadratic PC spectrum inferred for M + , we generate 50000 sample PC spectra {gi }50000 i=1 , where each spectrum is given by gi = {g0,i , . . . , g5,i }. Secondly, using the joint posterior, π(f0 , . . . , fPF ), of the PC spectrum inferred for F of a target transport coefficient, we generate 50000 corresponding sample PC spectra {fi }50000 of the appropriate order (recall, in fact, that a quadratic model is used for F ++ and i=1 +− µ F , a linear model for F , and a cubic for F −− ). For any given sample spectrum, gi and fi , we can compute one corresponding sample correlation coefficient Θi (M + , F ). For any given transport coefficient of interest, F µ , F ++ , F −− , and F +− , we can repeat this procedure for all the corresponding 50000 sample PC spectra. This yields a set of 50000 sample correlation coefficients for each combination. Figure 12 shows the histograms of the 50000 samples obtained for Θ(M + , F µ ) (a), Θ(M + , F ++ ) (b), Θ(M + , F −− ) (c), and Θ(M + , F +− ) (d). Panel (a) shows that Θ(M + , F µ ) is mainly negative. This result is in fact quite intuitive, as one expects the ionic flux to decrease when the viscosity increases. Panel (b) shows that Θ(M + , F ++ ) is tightly
18 −4
10
x 10
0.012 0.01 F ++ [mm 2 /s]
F µ [Pa · s]
9.8 9.6 9.4 9.2
0.006 0.004 0.002
9 0
25 M
+
50 [nm/ns]
75
0 0
100
0.012
0.012
0.01
0.01
0.008 0.006 0.004 0.002 0 0
25
(a)
F +− [mm 2 /s]
F −− [mm 2 /s]
0.008
50 M + [nm/ns]
75
50 M + [nm/ns]
75
100
(b)
0.008 0.006 0.004 0.002
25 M
+
50 [nm/ns]
75
0 0
100
(c)
25
100
(d)
FIG. 13. Scatter plots of the 10000 model predictions obtained using the MAP estimates of the PC spectra for (M + , F µ ) (a), (M + , F ++ ) (b), (M + , F −− ) (c), and (M + , F +− ) (d).
peaked around ∼ 0.96, and this is also anticipated because in the present regime one expects the flux of Na+ to be strongly correlated with its diffusivity. Panels (c,d) reveal that the correlations between M + and F −− , and between M + and the cross-diffusivity, F +− , weak and mostly positive. On the one hand, these results confirm the intuitive observation that the Na+ -conductance (or flux) is minimally affected by the diffusivity of the negative ions, F −− . On the other hand, however, we expected a stronger correlation between M + and the cross diffusivity, F +− . The weak correlation revealed here is probably a result of the dominant correlation found with F ++ shown in panel (b). Also, in the dilute limit, we would expect the correlation between M + and F −− to be practically zero, since, in this case, the N a+ ion flux would be independent of D−− . To further characterize the correlations, we perform the following sampling analysis. First, we generate 10000 samples, {ξ i }10000 i=1 , of the uniform random variables, ξ = {ξ1 , ξ2 }, defining the stochastic space. Secondly, these samples are used to compute the corresponding set of model predictions, {Mi+ }10000 i=1 , by using the MAP estimates of the PC coefficients, whose values can be obtained from the corresponding posterior density π(g0 , . . . , g5 ). We repeat this sampling procedure for the PCe of each transport coefficient, F , using the corresponding MAP estimates PC 10000 spectra, yielding the following four sets of model predictions Fiµ , Fi++ , Fi−− , Fi+− i=1 . Figure 13 shows the scatter plots of the samples obtained for (M + , F µ ) (a), (M + , F ++ ) (b), (M + , F −− ) (c), and (M + , F +− ) (d). The plots are consistent with the observations made before using the histograms of the correlation coefficients. Specifically, Figure 13 (a) shows that M + is overall negatively correlated with the viscosity. Panel (b) shows a clear positive correlation between M + and F ++ , while panels (c,d) reveal a weak correlation. Note that these plots reflect a nonlinear mapping of the stochastic space (ξ1 , ξ2 ) to the spaces (M + , F µ ), (M + , F ++ ), (M + , F −− ) and (M + , F +− ), namely due to the fact that transformation applied by M + is quadratic.
VI.
CONCLUSIONS
This article focused on quantifying parametric uncertainty stemming from a subset of potential parameters in MD simulations of concentration driven ionic flow through a silica nanopore. The model involved a silica nanopore connecting two reservoirs containing a solution of sodium and chloride ions in water, and an ad hoc concentration
19 control algorithm was implemented to create a counter flow of ions through the nanopore. Parametric uncertainty was introduced in two potential parameters, namely the Lennard-Jones energy parameters of both ions, εN a+ and εCl− . The uncertain ranges adopted for εN a+ and εCl− were based on current state of knowledge about these parameters, which reflects a wider uncertainty range for εN a+ than for εCl− . The heterogeneous nature of the system, combined with the impact of inherent thermal fluctuations, represented a key complexity of this study. PC representations built using Bayesian regression were exploited to suitably isolate the effect of the intrinsic noise, and to capture the effect of parametric uncertainty on the conductance of both ions. The results revealed weak correlations in the PC spectra obtained for both ionic conductances, and shapes closely resembling Gaussian, suggesting that the posterior distributions can also be suitably characterized in terms of marginalized densities. A Bayes factor analysis was applied, indicating that a second-order PC regression function is the best model to represent the Na+ conductance data, while a constant-value is the best model for the Cl− data. Due to the limited number of replicas considered, the trend of GCl− with respect to εCl− was obscured by the substantial noise level. On the other hand, the trend of GN a+ with respect to εN a+ was evident because the range of variation of εN a+ was sufficiently large for the inference analysis to overcome the impact of inherent fluctuations. Finally, the analysis explored the correlations between the ionic conductances and bulk transport coefficients, namely the viscosity and the collective diffusivities, which provided additional insight into the UQ results and the observed trends. Together with the experiences in Part I, this work provided detailed demonstration of the efficacy of the adopted UQ methods in characterizing parametric uncertainty and intrinsic noise in complex MD simulations. Efforts are underway to extend the present work to tackle broader problems, including the analysis of higher dimensional random spaces, as well as coupling stochastic MD predictions to mesoscale and continuum models. ACKNOWLEDGMENTS
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. Appendix A: Green-Kubo estimates of bulk transport coefficients
The Green-Kubo formulae [40–45] allow for the calculation of bulk transport coefficients from equilibrium molecular dynamics simulations. This linear response theory relates the integral of the correlation of the appropriate flux to an estimate of a given transport coefficient. The fluxes are phase variables, i.e. functions of the atomic positions xi and velocities vi . For instance, the dynamic viscosity, µ, can be estimated as Z ∞D E V µ = ς(0) · ς(t) dt (A1) 3kB T 0 where ς(t) is the deviatoric stress at time t, in a homogeneous system with temperature T , volume V , and with kB denoting the Boltzmann’s constant. The deviatoric stress {ς1 , ς2 , ς3 } = {σ12 , σ13 , σ23 } is zero in equilibrium and can be calculated from the virial. The virial expression for the stress σ [46–49] is : 1 X 1 X mi v i ⊗ v i + fij ⊗ xij σ= V 2! ij i 1 X + fijk ⊗ (xij + xik ) + . . . , (A2) 3! ijk
∂ φ2 (xi , xj ) are where xij ≡ xi − xj denotes the atomic position difference between the i-th and j-th atoms, fij ≡ − ∂x i the forces related to the pair potential φ2 , and
fijk ≡ −
∂ φ3 (xi , xj , xk ) ∂xi
20 are the forces related to the triplet potential φ3 , etc. Likewise, the (mutual) diffusion coefficients can be estimated via [50, 51]: Z N(a) N(b) ∞ ab hh(a) (0) · h(b) (t)i dt, D = N 0
(A3)
where N(a) is number of species a and N is the total number, and the dynamic quantity h(a) (t) of species a, is computed with respect to the solvent (species 0 is H2 O)) as h(a) =
X 1 1 X vα − vα , N(a) N(0) α∈G(a)
(A4)
α∈G(0)
where vα is the atom velocity, G(a) is the group containing all atoms of species (a) and N(a) is its size. Since there are many versions of these coefficients related to Fick’s law or the Maxwell-Stefan equations, we chose the simplest to compute since only the correlation with the pore flow is of interest. Note the correlation length of the velocity-velocity correlations used in this study is 16 ps, which is considerably longer than in Ref. [51].
[1] D. Frenkel and B. Smit, Understanding Molecular Simulation, Second Edition: From Algorithms to Applications, 2nd ed. (Academic Press, 2001). [2] A. Brodsky, Chemical Physics Letters 261, 563 (1996). [3] B. Guillot, J. Mol. Liq. 101, 219 (2002). [4] K. A. Dill, T. M. Truskett, V. Vlachy, and B. Hribar-Lee, Annual Review of Biophysics and Biomolecular Structure 34, 173 (2005). [5] A. Wallqvist and R. D. Mountain, “Molecular models of water: Derivation and description,” in Reviews in Computational Chemistry (John Wiley & Sons, Inc., 2007) pp. 183–247. [6] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, and J. Hermans, Intermolecular Forces , 331 (1981). [7] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 (1983). [8] F. Rizzi, R. Jones, B. Debusschere, and O. Knio., Physics of Fluids submitted (2012). [9] F. Rizzi, H. Najm, B. Debusschere, K. Sargsyan, M. Salloum, H. Adalsteinsson, and O. Knio., Multiscale Model. & Simul. in press (2012). [10] F. Rizzi, H. Najm, B. Debusschere, K. Sargsyan, M. Salloum, H. Adalsteinsson, and O. Knio., Multiscale Model. & Simul. in press (2012). [11] N. Wiener, Am. J. Math. 60, 897 (1938). [12] R. Ghanem and P. Spanos, Stochastic Finite Elements: A Spectral Approach (Springer Verlag, New York, 1991). [13] O. P. Le Maˆıtre and O. M. Knio, Spectral Methods for Uncertainty Quantification (Springer, New York, 2010). [14] R. Ghanem and A. Doostan, J. Comput. Phys. 217, 63 (2006). [15] K. Sargsyan, B. Debusschere, H. N. Najm, and Y. Marzouk, J. Comput. Theor. Nanosci. 6, 2283 (2009). [16] Note that the other non-bonded contribution originates from Coulombic forces between charged atoms. [17] M. W. Mahoney and W. L. Jorgensen, J. Chem. Phys. 112, 8910 (2000). [18] P. E. M. Lopes, V. Murashov, M. Tazi, E. Demchuk, and A. D. MacKerell, The Journal of Physical Chemistry B 110, 2782 (2006). [19] E. R. Cruz-Chu, A. Aksimentiev, and K. Schulten, The Journal of Physical Chemistry B 110, 21497 (2006). [20] C. D. Lorenz, P. S. Crozier, J. A. Anderson, and A. Travesset, The Journal of Physical Chemistry C 112, 10222 (2008), http://pubs.acs.org/doi/pdf/10.1021/jp711510k. [21] M. Patra and M. Karttunen, ArXiv Physics e-prints (2002), arXiv:physics/0211059. [22] W. L. Jorgensen and J. D. Madura, Mol. Phys. 56, 1381 (1985). [23] R. Ghanem, Comput. Methods Appl. Mech. Engrg. 158, 199 (1998). [24] O. P. Le Maˆıtre, M. T. Reagan, H. N. Najm, R. G. Ghanem, and O. M. Knio, J. Comput. Phys. 181, 9 (2002). [25] R. Ghanem, Comput. Methods Appl. Mech. Engrg. 168, 19 (1999). [26] R. Ghanem, J. Eng. Mech. 125, 26 (1999). [27] O. P. Le Maˆıtre, O. M. Knio, H. N. Najm, and R. G. Ghanem, J. Comput. Phys 173, 481 (2001). [28] M. Reagan, H. Najm, R. Ghanem, and O. Knio, Combustion and Flame 132, 545 (2003). [29] L. Fej´er, Bull. Amer. Math. Soc. 39, 521 (1933). [30] Multidimensional Fej´er grids are contained within the hypercube (−1, 1)d , where d is the dimension of the random vector ξ, and can be readily obtained by tensorization of the one-dimensional node points. In one dimension, the number of Fej´er points is given by nℓ = 2ℓ − 1, where ℓ is a desired level of resolution, and the points are obtained from the abscissae of the maxima of the Chebychev polynomials. The Fej´er grids are nested, i.e. the points at a resolution level, ℓ, also appear at higher levels ℓ′ > ℓ.
21 [31] A Jeffreys prior is a non-informative prior distribution that is proportional to the square root of the determinant of the Fisher information, and has the key feature that it is invariant under reparameterizations. This makes it particularly suitable for use with scale parameters. [32] H. Haario, E. Saksman, and J. Tamminen, Bernoulli 7, 223 (2001). [33] Y. F. Atchad´e and J. S. Rosenthal, Bernoulli 11, 815 (2005). [34] C. Andrieu and E. Moulines, Ann. Appl. Prob. 16, 1462 (2006). [35] G. O. Roberts and J. S. Rosenthal, J. Appl. Prob. 44, 458 (2007). [36] G. O. Roberts and J. S. Rosenthal, J. Comput. and Graph. Stat. 18, 349 (2009). [37] R. E. Kass and A. E. Raftery, Journal of the American Statistical Association 90, 773 (1995). [38] Note that the “collective” and self-diffusion are not the same quantity. Self-diffusion is a measure of the intrinsic mobility of a particle, and, in fact, is of interest for single component systems and Brownian motion. On the contrary, “collective” diffusion is relevant in multicomponent systems (e.g. particles in a solvent, like in the present case), since it accounts for the coupling between the motion of different particles, and reflects the contribution stemming from dynamic cross-correlations. To account for these important phenomena, a multicomponent (or collective) description of the diffusion process is required, which exploits linear combinations of the gradients of all concentrations (or chemical potentials) involved. This requires the knowledge of a full matrix of binary diffusion coefficients. In the present case, the ion-ion and ion-water interactions play a key role and, thus, this “collective” diffusion approach is necessary and appropriate. [39] The covariance, Cov(X, Y ), of two random variables, X and Y , is given by the expression appearing at the numerator in equation 10. The covariance provides a measure of the relative change between the two random variables, X, Y . If an increasing trend of X is associated with an increasing trend in Y , the covariance is a positive number. Conversely, when an increase in X is associated with a decrease in Y , the covariance is negative. Hence, the sign of the covariance shows the tendency in the linear relationship between the variables. [40] L. Onsager, Phys. Rev. 37, 405 (1931). [41] L. Onsager, Phys. Rev. 38, 2265 (1931). [42] M. S. Green, J. Chem. Phys. 22, 398 (1954). [43] R. Kubo, J. Phys. Soc. Jpn. 12, 570 (1957). [44] R. Kubo, M. Yokota, and S. Nakajima, J. Phys. Soc. Jpn. 12, 1203 (1957). [45] R. Zwanzig, J. Chem. Phys. 40, 2527 (1964). [46] R. J. E. Clausius, Philosophical Magazine 40, 122 (1870). [47] J. C. Maxwell, Transactions of the Royal Society Edinborough XXVI, 1 (1870). [48] J. H. Irving and J. G. Kirkwood, J. Chem. Phys. 18, 817 (1950). [49] K. K. Mandadapu, R. E. Jones, and P. Papadopoulos, Phys. Rev. E 80, 047702 (2009). [50] Y. Zhou and G. H. Miller, The Journal of Physical Chemistry 100, 5516 (1996). [51] D. R. Wheeler and J. Newman, The Journal of Physical Chemistry B 108, 18353 (2004).
