I Introduction
Molecular dynamics (MD) simulations are powerful computational methods for investigating the microscopic origins of the macroscopic behavior of soft materials such as viral capsids, lipid bilayers, confined electrolytes, selfassembled colloids, and polymers [1, 2, 3, 4]. These simulations furnish molecularlevel mechanisms that drive structure and property control in soft materials while isolating interesting regions of the material design space to aid experimental exploration and discovery. Despite the use of parallel computing techniques, simulations of soft materials often incur high computational costs. Recent years have seen a surge in the integration of machine learning (ML) methods with simulations to reduce their computational costs, enhance their predictive power, and expedite the analysis of highdimensional output data [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]. ML has been used to develop efficient force fields [10, 11, 12], reduce highdimensional simulation data to isolate molecularlevel mechanisms [7, 13], and develop surrogates to accurately predict simulation outcomes and expedite the exploration of the material design space [14, 15, 16, 17]. The last set of applications, where surrogates are designed for learning the relationship between the input variables and simulation outputs, is the subject of this paper.
Machine learning has been used to design surrogates to predict outputs of a wide variety of simulations [14, 15, 16, 17, 20, 21]. Deep neural networks trained on simulation data accurately predicted adsorption equilibria for different thermodynamic conditions [16]. ML surrogate for ab initio MD simulations bypassed the explicit time evolution of the particle trajectories to predict the dissociation timescale of compounds [15]
. Convolutional neural network based ML “emulators” have been developed to predict outputs of simulations in different scientific domains including climate science and biogeochemistry
[17]. Deep neural networkbased denoising autoencoder was trained to successfully predict the temporallyaveraged radial distribution function of LennardJones fluids from a single snapshot of the radial distribution function generated in MD simulations
[14]. A convolutional neural networkbased surrogate model was trained on data produced by MD simulation of polymer materials to predict filled rubber material properties that take the kinetic information of filler morphology as input [22].We recently introduced ML surrogates for MD simulations of soft materials [23, 24, 25, 26]. We demonstrated that artificial neural network based regression models trained on data from completed MD simulations of softmatter systems can successfully predict the relationships between the input parameters and the simulation outcomes, bypassing part or all of the explicit evolution of the simulated components. The approach was illustrated using simulations of electrolyte solutions confined by material surfaces. A key goal of the computational studies of confined electrolytes is to establish the links between the output distributions of electrolyte ions and the input electrolyte attributes for diverse solution conditions. These links provide a reliable guide to the regions of the material design space where significant changes in the structural organization of ions are expected, which can aid the experimental exploration and design of electrolytebased materials [27, 28, 29, 30, 31, 32]. The ML surrogate was trained to learn the relationship between the output distribution of positive ions and input variables characterizing the electrolytes containing positive and negative ions of the same size confined by uncharged surfaces. The ionic density profiles predicted by the surrogate were found to be in excellent agreement with MD simulations [23, 24]. Further, the inference time associated with the predictions was less than the corresponding MD simulation runtime.
In this paper, we extend the design and application of ML surrogates for MD simulations of soft materials in many significant ways. We show that statistical uncertainties present in the outputs of MD simulations can be utilized to train deep neural networks and design ML surrogates with higher accuracy and generalizability. We design soft labels for the simulation outputs by incorporating the uncertainties in the estimated average output quantities, and introduce a modified loss function that leverages these soft labels during training to significantly reduce the surrogate prediction error for input systems in the unseen test data. Using such soft labels for the ground truth data bypasses the need to explicitly expand the training dataset size in order to generalize the surrogate to learn complex relationships between input soft material attributes and output structural quantities.
The approach is illustrated with the design of a surrogate for MD simulations of confined electrolytes that predicts the distribution of positive and negative ions in confinement. The complexity of the relationship between input electrolyte attributes and output ionic structure is significantly enhanced by including ionic size asymmetry, charged surfaces, and a broad range of concentrations in the electrolyte model [33] (the surrogate developed in our previous paper [24] only predicted the distribution of positive ions for electrolytes with sizesymmetric ions confined by uncharged interfaces). We show that the surrogate, trained using the modified loss function that leverages the statistical uncertainties in the output ionic distributions, predicts the number densities of positive and negative ions in excellent agreement with the ground truth MD simulation results.
The high accuracy and small inference times (a factor of smaller than the corresponding simulation time) associated with the surrogate predictions of number density profiles enable quick access to structural quantities such as the net charge density and the integrated charge [34, 35]. The surrogateenabled generation of output ionic structure for millions of input electrolyte systems within minutes facilitates the implementation of complex analysis tasks such as sensitivity analysis that help distill the contributions of the input ionic attributes to the output ionic structure.
The rest of this paper is organized as follows. Section II describes the details of designing the ML surrogate. Section III discusses the accuracy characteristics and generalizability of the ML surrogate and presents the results of its application in confined electrolytes. Section IV presents conclusions.
Ii Machine Learning Surrogates for Molecular Dynamics Simulations
Figure 1 shows a system overview of the approach to design and use ML surrogates for MD simulations of soft materials [23, 24]. We first run MD simulations for different model variables (input) characterizing the softmatter system, and save the desired simulation outcomes (output) for training the ML surrogate, which occurs after a set number of successful simulation runs. Once trained, the ML surrogate is used to predict the relationship between the input variables and the output structural properties. Figure 1 highlights the modified loss function introduced in this paper to train the surrogate by leveraging the statistical uncertainties present in the MD simulation output data.
In the specific example of designing the surrogate for MD simulations of confined electrolytes, the input variables are confinement length, electrolyte concentration, ionic sizes, and surface charge, and the structural properties are the ionic density profiles. We begin by briefly describing these input variables and output quantities. A more detailed description of the model variables and the simulation methods and outputs can be found in our previous papers [33, 36].
Iia Input Variables and Output Quantities
We consider a monovalent electrolyte in water confined by two planar interfaces at room temperature K. The planes at and represent the location of the interfaces with being the interfacial separation ( corresponds to the midpoint between the interfaces). Each interface is characterized with a uniform charge density . A coarsegrained model [37, 38, 39, 40, 36] is employed to describe the electrolyte solution. Water is modeled as an implicit solvent with a relative dielectric permittivity of . The positivelycharged ions (cations) and negativelycharged ions (anions) associated with the monovalent electrolyte are modeled as spheres with hydrated sizes and respectively. An appropriate number of counterions, modeled as cations of the same diameter and charge, are included in the confinement to ensure electroneutrality.
Five input parameters: , electrolyte concentration , , , and , are used to design and train the ML surrogate (Table I). is defined as , where is the number of anions and is the volume of the simulation box. We note that not all ionic attributes and solutions conditions that are expected to alter the output ionic structure are considered as tunable input variables. For example, ion valencies (set to ), temperature (set at 298 K) and solvent permittivity (set to 80) are held fixed across all the simulations.
All MD simulations are performed using LAMMPS [41] in an NVT ensemble at K. Ionion and ioninterface steric interactions are modeled using LennardJones potentials [36], and all electrostatic interactions are modeled using Coulomb potentials whose longrange is properly treated with Ewald sums [42]. Each system is simulated for 1 ns to reach equilibrium with a timestep of femtosecond. After equilibration, systems are simulated for at least ns, and trajectory data are collected every ps.
Samples from the trajectory data are used to compute the average number densities and of cations and anions, respectively, together with the associated statistical uncertainties and . Note that due to the planar symmetry and the homogeneouslycharged interfaces, the ionic distributions vary only in the direction perpendicular to the interfaces, and hence are functions of a single variable . of cations and of anions form the output of the ML surrogate. These number densities are used to derive quantities such as the net charge density :
(1) 
and the integrated charge :
(2) 
Note that is defined using the left planar interface as the reference surface and is meaningful for , i.e., for up to the center of the confined region [33].
Input Variable  Range 

Interfacial separation ()  4 – 5 nm 
Electrolyte concentration ()  0.1 – 2.0 M 
Cation diameter ()  0.2 – 0.63 nm 
Anion diameter ()  0.2 – 0.63 nm 
Surface charge density ()  0.01, 0.015, 0.02 
IiB Data Generation and Preprocessing
Table I shows the ranges over which the input parameters , , , , and are tuned to generate the dataset for training and testing the ML surrogate. By sweeping over a few discrete values for each input variable within the associated range, 4,050 unique input configurations are created. For each of these configurations, MD simulations are run and the converged distributions of positive and negative ions are extracted as output using the ion trajectory samples collected during the MD simulation. These distributions are specified by computing the average population densities of ions at points (bins) within the confinement region (length) . The trajectory samples are also used to compute the associated statistical uncertainties in the extracted average number densities of cations and anions. Each MD simulation is performed for over nanoseconds and takes hours using LAMMPS with 96 cores. Generating the entire dataset (size GB) took approximately 20 days, including the queue wait times on the Indiana University BigRed3 supercomputing cluster.
The ML surrogate is trained to predict the average number densities of cations and anions for all bins. Thus, the surrogate makes a total of predictions to quantify the output ionic distributions. The entire dataset is separated into training and testing sets using a ratio of , resulting in training and testing datasets of 3240 and 810 configurations, respectively. A minmax normalization filter is applied to normalize the input data at the preprocessing stage.
IiC Feature Extraction and Regression
Figure 2 shows the artificial neural network (ANN) architecture employed in the ML surrogate to implement the regression and prediction of the desired output quantities. The ANN has 2 hidden layers, similar to the surrogate designed in our earlier work [24]
. The regression process determines the weights and biases in these hidden layers following an error backpropagation algorithm, implemented via a stochastic gradient descent procedure. This process uses a training dataset and an appropriate loss function for error computation and backpropagation to update the weights and biases after each batch of input data is regressed through the network by a simple forward prediction.
A common choice for the loss function in regression applications is the mean square error (MSE) between the ground truth data and the predictions :
(3) 
where denotes the input configuration and is the batch size. Note that and
are vectors associated with input system
whose size is equal to the number of output dimensions. In our earlier paper, we employed as the loss function [24] to train the ANN model to learn the density profiles of positive ions for electrolyte systems involving sizesymmetric ions and uncharged interfaces. In the present study, the ANN model is trained to learn a more complex inputoutput relationship due to the consideration of predicting both cation and anion density profiles for electrolytes involving sizeasymmetric ions, charged interfaces, and a wider range of concentrations. This rise in complexity is reflected in the larger errors in surrogate predictions for the test data when the ANN model is trained using (see Section IIIA). These errors can be reduced by increasing the number of training data samples, which requires running more MD simulations. Here, we adopt a different approach that bypasses the need to run these additional simulations. We develop a modified loss function that leverages the statistical uncertainties in the ground truth data to reduce the prediction errors and increase the generalizability of the ANN model.IiC1 Modified Loss Function
Many simulation outcomes such as density profiles of ions are associated with statistical uncertainties. These uncertainties are typically extracted as the standard error (standard deviation from the mean) using the samples (e.g., ion positions) collected during the simulation. For example, for a given input configuration, we extract the average (mean) number density of cations,
, and the associated statistical uncertainties , where the latter are given as:(4) 
Here, is the cation distribution computed using the trajectory sample , and is the number of samples. Trajectory samples are a set of ion positions collected every simulation time steps. Generally, includes errors resulting from the time discretization and the errors associated with the sampling process. The standard error captures the expected variations in the mean . By only considering the mean as the ground truth data while training a deep learning model, standard loss functions such as do not leverage the information encoded in regarding the variations in the ground truth data.
Our key idea is to utilize the uncertainties associated with the ground truth data to “informate” the training of the ANN model in order to enhance its generalizability when tested on unseen input samples. To implement the idea, we modify the standard MSE loss function by replacing the ground truth data vector
with an associated set of normal distributions having mean vector
and variance vector
, where is the uncertainty vector corresponding to . In other words, we replace the “hard” labels with soft labels , which are randomly sampled from the aforementioned normal distributions. For example, for each value of , we replace with samples drawn from a normal distribution having mean and variance . In general, the modified loss function can be written as:(5) 
where
is a random variable sampled from the normal distribution
. The backpropagated error through the ANN model obtained as the partial derivative of with respect to the model output is(6) 
Thus, enables the consideration of standard errors associated with the mean ground truth data to update the weights and biases when backpropagating through the network.
The use of soft labels for describing the ground truth data facilitates a sampling mechanism driven by the normal distributions that implicitly expands the dataset with more samples as the model undergoes training over many epochs. In other words, we enable the proliferation of samples without explicitly expanding the training dataset. Here onward, we refer to the model trained with the modified loss function
as ModNN and a baseline model trained with as BaseNN.We note that research in the applications of ML in computer vision and speech recognition has shown that uncertainties in the ground truth labels can be harnessed toward designing improved ML models, in particular, for classification tasks
[43, 44, 45, 46, 47]. For example, ground truth labels were represented as soft labels using a probability distribution to train neural networks
[43][45]for image classification tasks. Also, support vector machines for classifying voice quality samples yielded higher accuracy when trained on
fuzzy output (soft labels) compared to standard output [47].Aung et. al. studied the comparison between soft labels and hard labels for training recognition models as classifiers and regressors for computer vision applications [44]. Their findings showed that soft labels can improve the accuracy for both classification and regression tasks due to the enhanced regularization effect. Other studies have explored the connection between improved regularization capabilities via the use of soft labels and dropout regularization [48]. The impact of noise in input or output data on the generalization of ML models has also been studied [49, 50]. Our approach introduces an alternate method based on the modification of the loss function to utilize soft labels during the training of deep neural networks for regression tasks.
IiC2 Surrogate Training
The implementation, training, and testing of the artificial neural network (ANN) model were programmed using scikitlearn, Keras, and TensorFlow ML libraries
[51, 52, 53]. Scikitlearn was used for grid search and feature scaling, Keras was used to save and load models, and TensorFlow was used to define and train the neural network. Training the ANN model involves an appropriate selection of hyperparameters such as the number of first hidden layer units
, the number of second hidden layer units , batch size , and the number of epochs . is a hyperparameter of the stochastic gradient descent algorithm that controls the number of training samples allowed to pass through the ANN before updating its trainable parameters (e.g., weights, biases). is a hyperparameter that controls the number of complete passes made through the entire training dataset.By performing a grid search, we observed that the hyperparameters are optimized to , , , and for both BaseNN and ModNN models. The Adam optimizer was used to optimize the error backpropagation process. Using a trialanderror process, the learning rate of the Adam optimizer and the dropout rate in the dropout layer were set to and , respectively. The weights in the hidden layers and the output layer were initialized for better convergence using a Xavier normal distribution at the beginning. This distribution is characterized with 0 mean and variance, where and are the input and output sizes of the hidden layers, respectively [54]. Prototype implementation written using Python/C++ is available on GitHub [55], and a software application utilizing the surrogate is deployed in nanoHUB [56].
Iii Results and Discussion
Iiia Surrogate Accuracy: BaseNN vs ModNN Models
We begin by comparing the generalizability of the ML surrogate trained using the modified loss function (ModNN model) vs the surrogate trained using the standard loss function (BaseNN model). The generalization ability can be discussed in terms of overfitting characteristics and accuracy assessed by examining the error incurred on surrogate predictions for inputs in the unseen test data.
We first compute the test loss for each epoch as the average mean square error (MSE) incurred in making predictions describing the cation and anion density profiles associated with input systems in the test data:
(7) 
Here is the prediction vector for the input system comprising elements, and is the corresponding ground truth vector. Note that is the prediction of the ion number density for the input system .
Figure 3 compares the test loss (MSE) for the ModNN vs the BaseNN models. MSE of ModNN model is smaller than the BaseNN model for almost all epochs. The test loss decreases to an optimal value of near 10000 epochs for the BaseNN model, but goes down further to within 10000 epochs for the ModNN model. Within the same number of epochs, we find the optimal training loss for the BaseNN and ModNN models to be and respectively. Similar levels of test and train loss values suggest that both models are not overfitted near 10000 epochs [57]. The test loss for the BaseNN model increases gradually when the model is trained beyond 10000 epochs, signaling the crossover to the overfitting territory. In contrast, the test loss for the ModNN model does not cross over to overfitting territory even after 20000 epochs of training.
Figure 3 (inset) breaks down the test loss into contributions emerging from predictions near the interfaces (within nm from either interface) and the predictions in the bulk near the center of the confinement ( nm). The optimal test loss for predictions near the interfaces is higher than that for predictions in the bulk for both BaseNN and ModNN models. This difference stems from the relatively higher complexity of the relationship between the interfacial ionic structure and the electrolyte attributes.
MSE for BaseNN model predictions near the interfaces and in the bulk decreases initially with increasing number of epochs, reaching lowest values of (interfaces) and (bulk) around 10000 epochs, and then increases gradually on further training. MSE for ModNN model predictions near the interfaces and in the bulk decreases with increasing number of epochs, reaching (interfaces) and (bulk) around 10000 epochs, and does not increase on further training. Compared to the BaseNN model, the ModNN model reduces the prediction error near the interfaces and in bulk by over and over respectively.
We now examine the accuracy associated with the optimal BaseNN and ModNN models produced by selecting weights and biases checkpointed at 9670 epochs using early stopping criteria. The accuracy is measured for each of the ion density predictions associated with an input system by computing the root mean square errors (RMSE). For the prediction, is the error associated with the prediction of the ion density at the bin:
(8) 
where, is the prediction (inference) of the ion number density at the bin for the input system specified by the index , is the corresponding ground truth, and . is evaluated by averaging over the errors associated with the bin for all the samples in the test dataset. It is useful to note that the values are related to the MSE extracted using Equation 7 at 9670 epoch via the equality: .
For each bin, 2 inferences are made: one for the cation density and the other for the anion density, such that half of the predictions () describe the cation density profile, while the other half describes the anion density profile. We collect these 2 groups of inferences into 2 separate RMSE sets, and , with . The bin corresponds to the discretization of the confinement length and is correlated with the distance from the left interface. Thus, and also indicate the variation of the error in the prediction of cation and anion densities as a function of the location within the confinement.
Figure 4 shows and associated with the predictions made by the ModNN and BaseNN models for the cation and anion densities in the bin. For all , and values produced by the ModNN model are smaller than those produced by the BaseNN model. The average is 0.0015 and 0.00042 for the BaseNN and ModNN models respectively. The average is 0.0012 and 0.00041 for the BaseNN and ModNN models respectively. and are higher near the interfaces ( or ) compared to the bulk () for either models. This difference was also observed in the MSE values (Figure 3 inset), and can be attributed to the greater complexity in the variations of the ionic structure near the interfaces due to changes in electrolyte attributes.
The lower MSE values (Figure 3) for the ModNN model predictions for systems in the unseen test dataset show that the training of the surrogate using the modified loss function leads to a more optimal selection of the weights and biases compared to training the surrogate using the standard loss function . Further, the significantly smaller RMSE associated with the predictions made by the ModNN model (Figure 4) demonstrate that the surrogate exhibits higher generalizability when trained with . In the following sections, we use the ML surrogate trained using the modified loss function (ModNN model) to predict the ionic structure.
IiiB Ionic Density Profiles
We compare the cation and anion number density profiles predicted by the ML surrogate with the ground truth results obtained using MD simulations for input systems in the test dataset. Figure 5 shows the results for 4 input systems (I, II, III, IV) randomly selected from the test dataset. The systems are labeled by 5 input variables characterizing the ionic system (defined in Section IIA): , , , , . The 4 systems are: system I (4.8, 1.5, 0.415, 0.63, 0.02), system II (5, 2, 0.63, 0.63, 0.02), system III (4.4, 1.0, 0.415, 0.5225, 0.015), and system IV (4, 1.5, 0.5225, 0.415, 0.015).
For each system, the MLpredicted ion density profiles are in excellent agreement with the ground truth results extracted using MD simulations. A quantitative assessment of the accuracy of the predicted density profiles by the ML surrogate can be made by extracting the RMSE values for an input system for cation and anion density profiles:
(9) 
where, is the prediction of the cation or anion number density at the bin, is the corresponding ground truth, and is the total number of predictions per density profile. values associated with the cation density profiles for systems I, II, III, and IV are 0.0007, 0.001, 0.0003, and 0.0004 respectively. values associated with the anion density profiles for systems I, II, III, and IV are 0.0004, 0.0009, 0.0003, and 0.0005 respectively.
It is useful to provide additional metrics to assess the accuracy of the surrogate predictions for the density profiles. Following our previous paper [24], we compute the success rate associated with the input system as:
(10) 
where is the surrogate prediction of the ion density at the bin, and are respectively the corresponding ground truth result and error obtained using MD simulations, and is a step function defined as for ; for . Note that the sum in Equation 10 is over the entire set of predictions made using the ML surrogate, which includes predictions for both cation and anion densities. The success rate can be thought of as the number of instances when the ML prediction of density values for the given ionic system were found to be within the error bounds (Equation 4) produced by the associated MD simulation. For example, a success rate of implies that out of 1004 density predictions, 974 were within the associated uncertainties obtained from MD simulations.
The success rates for the prediction of density profiles for systems I, II, III, and IV are , , , and respectively. These values are similar to the success rates computed in our previous paper for ionic systems characterized with uncharged surfaces and sizesymmetric ions. Thus, the surrogate trained using the modified loss function yields high success rates despite the increase in the complexity of the ionic systems considered in this work. We note that the success rates using the BaseNN model (standard loss function ), are around for these 4 systems.
In addition to the high accuracy of the surrogate, the inference time associated with ML predictions of density profiles ( seconds) is a factor of smaller compared to the corresponding parallelized MD simulation runtime ( hours). The combination of high accuracy and small inference times associated with the surrogate predictions provides rapid access to quantities derived using the density profiles and facilitates the implementation of complex analysis tasks. We will discuss these results next.
IiiC Net Charge Density and Integrated Charge
The excellent agreement between the cation and anion number densities generated via ML surrogate and MD simulations paves the way for using the surrogate to derive the net charge density (Equation 1) and the integrated charge (Equation 2). The integrated charge provides a way to quantify the effects of ion accumulation and depletion near the surface on the screening of the surface charge. We define using the left planar interface as the reference surface, and plot it as a function of , where the latter denotes the distance from the left interface. Note that is meaningful for (or ), i.e., for a distance less than half the confinement length.
Figure 6 shows the integrated charge and the net charge density (inset) computed by the ML surrogate for the same aforementioned four systems (I, II, III, and IV). The corresponding ground truth results obtained using MD simulations are also shown. For each system, the integrated charge and the net charge density computed by the ML surrogate are in excellent agreement with the ground truth results.
In a recent study [33], using the structural properties and , we demonstrated the existence of two distinct regimes of screening behavior for different electrolyte systems as the electrolyte concentration increases from 0.1 M to 2.5 M. While a broad range of systems was generated and explored by tuning , , , and , the higher simulation costs limit the campaign to probe the material design space to hundreds of electrolyte systems. The availability of an accurate and fast ML surrogate enables a reliable access to variations in the ionic structure with continuous changes in the solution conditions for millions of electrolyte systems within minutes. This enables campaigns to probe the continuous input design space within the boundaries set by the initial simulation exploration used to design the surrogate.
Figure 7 shows an example of a surrogateled campaign generating the variation in the net charge density and the integrated charge with changing electrolyte concentration 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7, and 1.9 M for an electrolyte system with nm, nm, nm, and C/m. Results are shown for the left half of the confinement within 1.5 nm from the interface. Increasing leads to enhanced accumulation of cations and anions near the interface as evident by the rise in the magnitude of the peaks in near the interface. The pattern of the decay of to 0 as also changes as is increased. For low , exhibits a monotonic decay with a characteristic decay length decreasing sharply with increasing . However, for high , shows a nonmonotonic oscillatory decay with a decay length exhibiting a small rise with increasing . These findings produced by the ML surrogate track the behavior reported in our recent paper [33].
IiiD Sensitivity Analysis
The rich variations observed in the net charge density and the integrated charge as a function of the location within the confinement (Figure 6), and the dramatic changes in these quantities with increasing concentration (Figure 7) are linked to the complex organization of cations and anions within the confinement. For example, simulation snapshots extracted using the trajectory data indicate the emergence of structured layers of cations and anions near the interface when the concentration is increased to values M [33]. The high accuracy and small inference times associated with the ML surrogate predictions facilitate the implementation of complex analysis tasks for probing the observed modulations in the ionic structure and distilling the contributions of the attributes of cations and anions to the structural properties.
The surrogate can generate millions of cation and anion density profiles within minutes to yield smooth response surfaces describing variations of output structural properties to any changes in input electrolyte attributes. Leveraging this information, the surrogate can perform rapid sensitivity analysis and quantitatively assess the relative contributions of different ion attributes to the observed structural properties. Specifically, we demonstrate the use of the ML surrogate to perform sensitivity analysis and extract sensitivities of the response of the net charge density to variations in the cation size and the anion size for different locations within the confinement.
The sensitivity analysis can be performed using many methods [58, 59]. We use variancebased methods [60, 61, 59], wherein the variance of the output is decomposed into fractions that can be attributed to distinct input variables. We present the process for computing the sensitivity of the output quantity for the case of two input variables: and . The steps can be readily generalized to a large number of outputs and inputs. We begin by identifying the response surface connecting the distribution associated with the output to the distributions associated with the input variables , for each . Generally, are randomly sampled from distributions with mean values specified by a fixed pair of .
Let be the expectation value (mean) of . Let be the associated variance. In the variancebased methods, we compute the contribution of varying only one input variable (e.g., ) to the overall output variance , which can be described as an average over the variations in all the other input variables (e.g., ). We first define as the average of the output net charge density for a fixed sampled from a distribution (e.g., normal). This average is computed by summing over number of samples drawn from a distribution associated with the input variable . The process is repeated for all values, and for are generated, where denotes the total number of samples. The variance measuring the effect of varying only is the variance of this resulting output, :
(11) 
The sensitivity quantifying the contribution of varying only to the overall output variance is
(12) 
Following the same process and replacing with , the sensitivity quantifying the contribution of varying only to the overall output variance is:
(13) 
Note that .
We now extract the sensitivities and for the net charge densities representative of the distinct regimes of screening behavior found under low and high conditions. We select two electrolyte systems, one at low M and another at high M; other input variables for both systems are fixed at nm and . The mean values for the cation and anion diameter are taken to be nm and nm respectively. As the net charge density profiles for these two systems involve the assessment of for a large set of values, one can extract a total of pairs of sensitivities and . For the ease of illustration, we limit our analysis to 5 values within the left half of the confinement: (I) nm, (II) nm, (III) nm, (IV) nm, and (V) nm. These values are chosen as representative of the regions within confinement where ionic size effects are expected to influence distinct features of the ionic structure. Let denote the distance from the left interface. (I) and (II) represent the points of contact density of the cation () and anion () respectively, (III) represents a distance between and , (IV) represents a distance slightly greater than , and (V) represents a point near the center of the confinement, within the bulk region ( nm) for either systems.
We generate 1000 samples of cation sizes from the normal distribution nm, and 1000 samples of anion sizes from nm, thus creating a collection of 1 million electrolyte systems. Note that nm and nm. We then use the ML surrogate to predict the cation and anion number density profiles, and the associated net charge densities, for the entire set of 1 million samples. This process generates the response surfaces connecting the distribution associated with the output to the distributions associated with input variables for the aforementioned five values. Leveraging the response surfaces, we extract the sensitivities using Equations 12 and 13.
Figure 8 shows the response surfaces and associated sensitivities and for four values (I – IV) for the low M (Figure 8A) and high M (Figure 8C) electrolyte systems. We represent the sensitivities as percentages in the form of a pie chart next to the 3D plot of the response surface . The corresponding net charge density profiles for the low and high systems are shown in Figure 8B and Figure 8D respectively.
At the point of contact density of cations (case I: ), we find and for both low and high systems. This indicates that the sensitivity of to changes in ion size is almost entirely dominated by variations in the cation size at the distance away from the interface. This finding is expected because anions, which on average are larger compared to cations in the system considered, are excluded from the neighborhood of (Figure 8B and D), and therefore the steric effects influencing the ionic structure emerge entirely due to the cation size regardless of the electrolyte concentration.
The contribution of the size of the anions to the sensitivity of becomes prominent near the point of contact density of anions (case II: ). For M, we find , i.e., is much more sensitive to the variations in the anion size compared to those in the cation size. We attribute this to the depletion of anions resulting from the negativelycharged interface in the vicinity of , as evident by the associated anion density profile (Figure 8C). Small variations in the anion size alter the location of the anion contact density. Therefore anion density at changes with in a substantial way, making highly sensitive to variations in anion size compared to cation size. Results are different for M: is only marginally larger compared to . The anion contact density is not as sensitive to the anion size because of the strong steric ionion repulsion that pushes anions near the interface regardless of their size and repulsion from the negativelycharged interface.
At a point that represents a distance between and (case III: nm), we observe a clear switch in the sensitivity profile between the low and high systems. For the low c system, is much more sensitive to the changes in cation size compared to anion size (; ). Here, the accumulation of cations is higher than anions (Figure 8C). Changes in cation size lead to stronger variations in cation density compared to those in anion density resulting from changes in anion size. The higher sensitivity of to changes in anion size in the high system can be attributed to the vicinity of nm inhabiting the location of the anion peak density, which dominates the contribution to compared to the relatively depleted cation density (Figure 8D). These findings suggest that the peak of anion density is highly sensitive to the anion size.
For distance slightly greater than anion diameter from the interface (case IV: nm, nm), for the low system, indicating that the variations in cation and anion sizes contribute equally to the sensitivity of . This can be attributed to the weak modulations in the ionic structure at distances greater than the anion diameter, which can be identified as the onset of the bulk region for the low system (Figure 8C). Similar and values indicate that variations in cation and anion sizes contribute in equal amounts to changes in the bulk value of . For the same location within the confinement ( nm), for the high system, indicating that the modulations in are still strong at this distance from the interface. At nm, the high system has not reached the bulk region (also evident in Figure 8D). The switch in the sensitivity profile compared to case III, where and , originates from the presence of alternating layers of cations and anions near the interface. This structured region extends much farther into the center of the confinement compared to the low system. Only at distances greater than twice the anion diameter from the interface (case V: nm, nm), we find for the high system, signaling the onset of the bulk region.
Iv Conclusion
Output quantities produced by molecular dynamics simulations of soft materials are generally associated with statistical uncertainties. We showed that these uncertainties can be utilized to informate the training of deep neural networks for designing machine learning surrogates aimed at predicting the relationship between input variables and simulation outputs. The approach was illustrated with the design of a surrogate for molecular dynamics simulations of confined electrolytes to predict the complex relationship between the input electrolyte attributes and the output ionic structure. We demonstrated that the prediction error for samples in the unseen test data can be significantly reduced by utilizing a modified loss function that leverages the uncertainties in the output ionic distributions during training.
Using soft labels for the ground truth data facilitates a sampling mechanism that implicitly expands the dataset with more samples as the model undergoes training over many epochs, producing a surrogate that exhibits higher generalizability. Comparing the improvement in accuracy resulting from the use of the modified loss function with results obtained using surrogates trained by explicitly expanding the training dataset size will be a subject of future work. Further, we will explore the extension of these ideas to develop surrogates for simulations of other soft materials exhibiting complex inputoutput relationships.
The surrogate predictions for the cation and anion number density profiles were found to be in excellent agreement with the ground truth results produced using molecular dynamics simulations. The high accuracy and small inference times associated with the surrogate predictions enabled rapid derivation of the net charge density and the integrated charge using the ion number densities, and facilitated the implementation of complex analysis tasks such as sensitivity analysis. We demonstrated that the surrogate can rapidly extract sensitivities of the net charge density to cation and anion sizes, thus aiding in distilling the contributions of the attributes of cations and anions to the distinct ionic structure observed at low and high electrolyte concentrations.
Acknowledgments
This work is supported by the National Science Foundation through Awards 1720625 and DMR1753182, and by the Department of Energy through Award DESC0021418. Simulations were performed using the BigRed3 supercomputing system at Indiana University.
References
 [1] S. C. Glotzer, “Assembly engineering: Materials design for the 21st century (2013 pv danckwerts lecture),” Chemical Engineering Science, vol. 121, pp. 3–9, 2015.
 [2] J. Perilla, J. Hadden, B. Goh, C. Mayne, and K. Schulten, “Allatom molecular dynamics of virus capsids as drug targets,” The journal of physical chemistry letters, vol. 7, no. 10, pp. 1836–1844, 2016. [Online]. Available: https://pubs.acs.org/doi/10.1021/acs.jpclett.6b00517
 [3] N. E. Brunk and V. Jadhao, “Computational studies of shape control of charged deformable nanocontainers,” Journal of Materials Chemistry B, vol. 7, p. 6370, 2019. [Online]. Available: http://dx.doi.org/10.1039/C9TB01003C
 [4] V. Jadhao and M. O. Robbins, “Rheological properties of liquids under conditions of elastohydrodynamic lubrication,” Tribology Letters, vol. 67, no. 3, p. 66, 2019. [Online]. Available: https://doi.org/10.1007/s1124901911783

[5]
A. L. Ferguson, “Machine learning and data science in soft materials engineering,”
Journal of Physics: Condensed Matter, vol. 30, no. 4, p. 043002, 2017.  [6] M. Spellings and S. C. Glotzer, “Machine learning for crystal identification and discovery,” AIChE Journal, vol. 64, no. 6, pp. 2198–2206, 2018.
 [7] S. S. Schoenholz, E. D. Cubuk, E. Kaxiras, and A. J. Liu, “Relationship between local structure and relaxation in outofequilibrium glassy systems,” Proceedings of the National Academy of Sciences, vol. 114, no. 2, pp. 263–267, 2017. [Online]. Available: https://doi.org/10.1073/pnas.1610204114
 [8] A. Z. Guo, E. Sevgen, H. Sidky, J. K. Whitmer, J. A. Hubbell, and J. J. de Pablo, “Adaptive enhanced sampling by forcebiasing using neural networks,” The Journal of chemical physics, vol. 148, no. 13, p. 134108, 2018.
 [9] M. G. Wessels and A. Jayaraman, “Machine learning enhanced computational reverse engineering analysis for scattering experiments (crease) to determine structures in amphiphilic polymer solutions,” ACS Polymers Au, pp. 8581–8591, 2021.
 [10] J. Wang, S. Olsson, C. Wehmeyer, A. Perez, N. E. Charron, G. De Fabritiis, F. Noe, and C. Clementi, “Machine learning of coarsegrained molecular dynamics force fields,” ACS central science, 2019.
 [11] L. Casalino, A. C. Dommer, Z. Gaieb, E. P. Barros, T. Sztain, S.H. Ahn, A. Trifan, A. Brace, A. T. Bogetti, A. Clyde et al., “Aidriven multiscale simulations illuminate mechanisms of sarscov2 spike dynamics,” The International Journal of High Performance Computing Applications, p. 10943420211006452, 2021.
 [12] W. Jia, H. Wang, M. Chen, D. Lu, L. Lin, R. Car, W. E, and L. Zhang, “Pushing the limit of molecular dynamics with ab initio accuracy to 100 million atoms with machine learning,” in Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, ser. SC ’20, 2020. [Online]. Available: https://dl.acm.org/doi/abs/10.5555/3433701.3433707
 [13] J. Kadupitiya and V. Jadhao, “Probing the rheological properties of liquids under conditions of elastohydrodynamic lubrication using simulations and machine learning,” Tribology Letters, vol. 69, no. 3, pp. 1–19, 2021. [Online]. Available: https://doi.org/10.1007/s11249021014573

[14]
A. Moradzadeh and N. R. Aluru, “Molecular dynamics properties without the full trajectory: A denoising autoencoder network for properties of simple liquids,”
The journal of physical chemistry letters, vol. 10, no. 24, pp. 7568–7576, 2019. [Online]. Available: https://doi.org/10.1021/acs.jpclett.9b02820  [15] F. Häse, I. Fdez. Galván, A. AspuruGuzik, R. Lindh, and M. Vacher, “How machine learning can assist the interpretation of ab initio molecular dynamics simulations and conceptual understanding of chemistry,” Chem. Sci., vol. 10, pp. 2298–2307, 2019. [Online]. Available: http://dx.doi.org/10.1039/C8SC04516J
 [16] Y. Sun, R. F. DeJaco, and J. I. Siepmann, “Deep neural network learning of complex binary sorption equilibria from molecular simulation data,” Chemical science, vol. 10, no. 16, pp. 4377–4388, 2019.
 [17] M. Kasim, D. WatsonParris, L. Deaconu, S. Oliver, P. Hatfield, D. Froula, G. Gregori, M. Jarvis, S. Khatiwala, J. Korenaga et al., “Building high accuracy emulators for scientific simulations with deep neural architecture search,” arXiv eprints, pp. arXiv–2001, 2020. [Online]. Available: https://arxiv.org/abs/2001.08055
 [18] J. Kadupitiya, G. C. Fox, and V. Jadhao, “Machine learning for parameter autotuning in molecular dynamics simulations: Efficient dynamics of ions near polarizable nanoparticles,” The International Journal of High Performance Computing Applications, vol. 34, no. 3, pp. 357–374, 2020. [Online]. Available: https://journals.sagepub.com/doi/full/10.1177/1094342019899457
 [19] ——, “Deep learning based integrators for solving newton’s equations with large timesteps,” arXiv preprint arXiv:2004.06493, 2021. [Online]. Available: https://arxiv.org/abs/2004.06493
 [20] M. T. Degiacomi, “Coupling molecular dynamics and deep learning to mine protein conformational space,” Structure, vol. 27, no. 6, pp. 1034–1040, 2019.
 [21] A. Rahman, P. Deshpande, M. S. Radue, G. M. Odegard, S. Gowtham, S. Ghosh, and A. D. Spear, “A machine learning framework for predicting the shear strength of carbon nanotubepolymer interfaces based on molecular dynamics simulation data,” Composites Science and Technology, vol. 207, p. 108627, 2021.
 [22] T. Kojima, T. Washio, S. Hara, and M. Koishi, “Synthesis of computer simulation and machine learning for achieving the best material properties of filled rubber,” Scientific Reports, vol. 10, no. 1, pp. 1–11, 2020.
 [23] J. Kadupitiya, G. C. Fox, and V. Jadhao, “Machine learning for performance enhancement of molecular dynamics simulations,” in International Conference on Computational Science, 2019, pp. 116–130.
 [24] J. Kadupitiya, F. Sun, G. Fox, and V. Jadhao, “Machine learning surrogates for molecular dynamics simulations of soft materials,” Journal of Computational Science, p. 101107, 2020.
 [25] G. Fox, J. A. Glazier, J. Kadupitiya, V. Jadhao et al., “Learning everywhere: Pervasive machine learning for effective highperformance computation,” in IEEE IPDPS Workshops, 2019, pp. 422–429. [Online]. Available: https://doi.org/10.1109/IPDPSW.2019.00081
 [26] V. Jadhao and J. Kadupitiya, “Integrating machine learning with hpcdriven simulations for enhanced student learning,” in 2020 IEEE/ACM Workshop on Education for HighPerformance Computing (EduHPC). IEEE, 2020, pp. 25–34. [Online]. Available: https://ieeexplore.ieee.org/document/9308675
 [27] Y. He, D. Gillespie, D. Boda, I. Vlassiouk, R. S. Eisenberg, and Z. S. Siwy, “Tuning transport properties of nanofluidic devices with local charge inversion,” Journal of the American Chemical Society, vol. 131, no. 14, pp. 5194–5202, 2009.
 [28] S. J. Faucher, N. R. Aluru, M. Z. Bazant, D. Blankschtein, A. H. Brozena, J. Cumings, J. P. de Souza, M. Elimelech, R. Epsztein, J. T. Fourkas et al., “Critical knowledge gaps in mass transport through singledigit nanopores: A review and perspective,” The Journal of Physical Chemistry C, 2019.
 [29] H. B. Park, J. Kamcev, L. M. Robeson, M. Elimelech, and B. D. Freeman, “Maximizing the right stuff: The tradeoff between membrane permeability and selectivity,” Science, vol. 356, no. 6343, p. eaab0530, 2017.
 [30] J. R. Werber, C. O. Osuji, and M. Elimelech, “Materials for nextgeneration desalination and water purification membranes,” Nature Reviews Materials, vol. 1, no. 5, pp. 1–15, 2016.
 [31] Y. Levin, “Strange electrostatics in physics, chemistry, and biology,” Physica A: Statistical Mechanics and its Applications, vol. 352, no. 1, pp. 43 – 52, 2005.
 [32] J. W. Zwanikken and M. O. de la Cruz, “Correlated electrolyte solutions and ioninduced attractions between nanoparticles,” Phys. Rev. E, vol. 82, p. 050401, Nov 2010.
 [33] N. Anousheh, F. J. Solis, and V. Jadhao, “Ionic structure and decay length in highly concentrated confined electrolytes,” AIP Advances, vol. 10, no. 12, p. 125312, 2020. [Online]. Available: https://aip.scitation.org/doi/10.1063/5.0028003
 [34] R. Qiao and N. R. Aluru, “Charge inversion and flow reversal in a nanochannel electroosmotic flow,” Physical review letters, vol. 92, no. 19, p. 198301, 2004.

[35]
M. F. Dopke, J. Lutzenkirchen, O. A. Moultos, B. Siboulet, J.F. Dufrêche, J. T. Padding, and R. Hartkamp, “Preferential adsorption in mixed electrolytes confined by charged amorphous silica,”
The Journal of Physical Chemistry C, vol. 123, no. 27, pp. 16 711–16 720, 2019.  [36] Y. Jing, V. Jadhao, J. W. Zwanikken, and M. Olvera de la Cruz, “Ionic structure in liquids confined by dielectric interfaces,” The Journal of chemical physics, vol. 143, no. 19, p. 194508, 2015.
 [37] D. Boda, D. Gillespie, W. Nonner, D. Henderson, and B. Eisenberg, “Computing induced charges in inhomogeneous dielectric media: Application in a monte carlo simulation of complex ionic systems,” Phys. Rev. E, vol. 69, no. 4, p. 046702, Apr 2004.
 [38] R. Allen, J.P. Hansen, and S. Melchionna, “Electrostatic potential inside ionic solutions confined by dielectrics: a variational approach,” Phys. Chem. Chem. Phys., vol. 3, pp. 4177–4186, 2001.
 [39] S. Tyagi, M. Suzen, M. Sega, M. Barbosa, S. S. Kantorovich, and C. Holm, “An iterative, fast, linearscaling method for computing induced charges on arbitrary dielectric boundaries,” The Journal of Chemical Physics, vol. 132, no. 15, p. 154112, 2010.
 [40] K. Barros, D. Sinkovits, and E. Luijten, “Efficient and accurate simulation of dynamic dielectric objects,” The Journal of Chemical Physics, vol. 140, no. 6, p. 064903, 2014.
 [41] S. Plimpton, “Fast parallel algorithms for shortrange molecular dynamics,” Journal of Computational Physics, vol. 117, no. 1, pp. 1 – 19, 1995.
 [42] M. Deserno and C. Holm, “How to mesh up ewald sums. i. a theoretical and numerical comparison of various particle mesh routines,” The Journal of chemical physics, vol. 109, no. 18, pp. 7678–7693, 1998.
 [43] P. Smyth, M. C. Burl, U. M. Fayyad, and P. Perona, “Knowledge discovery in large image databases: Dealing with uncertainties in ground truth.” in KDD workshop, 1994, pp. 109–120.
 [44] A. M. Aung and J. Whitehill, “Harnessing label uncertainty to improve modeling: An application to student engagement recognition.” in FG, 2018, pp. 166–170.
 [45] S. Kumano, K. Otsuka, D. Mikami, M. Matsuda, and J. Yamato, “Analyzing interpersonal empathy via collective impressions,” IEEE Transactions on Affective Computing, vol. 6, no. 4, pp. 324–336, 2015.
 [46] H. Meng, A. Kleinsmith, and N. BianchiBerthouze, “Multiscore learning for affect recognition: the case of body postures,” in International Conference on Affective Computing and Intelligent Interaction. Springer, 2011, pp. 225–234.
 [47] S. Scherer, J. Kane, C. Gobl, and F. Schwenker, “Investigating fuzzyinput fuzzyoutput support vector machines for robust voice quality classification,” Computer Speech & Language, vol. 27, no. 1, pp. 263–287, 2013.
 [48] S. Wager, S. Wang, and P. S. Liang, “Dropout training as adaptive regularization,” Advances in neural information processing systems, vol. 26, pp. 351–359, 2013.

[49]
A. Krogh and J. A. Hertz, “Generalization in a linear perceptron in the presence of noise,”
Journal of Physics A: Mathematical and General, vol. 25, no. 5, p. 1135, 1992.  [50] K. Crammer, A. Kulesza, and M. Dredze, “Adaptive regularization of weight vectors,” Advances in neural information processing systems, vol. 22, 2009.
 [51] F. Chollet et al., “Keras,” 2015.
 [52] L. Buitinck et al., “Api design for machine learning software: experiences from the scikitlearn project,” arXiv:1309.0238, 2013.
 [53] M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard et al., “Tensorflow: a system for largescale machine learning.” in OSDI, vol. 16, 2016, pp. 265–283.

[54]
X. Glorot and Y. Bengio, “Understanding the difficulty of training deep
feedforward neural networks,” in
Proceedings of the thirteenth international conference on artificial intelligence and statistics
, 2010, pp. 249–256.  [55] GitHub, “Repository nanoconfinementmd in softmaterialslab,” 2021. [Online]. Available: https://github.com/softmaterialslab/nanoconfinement
 [56] K. Kadupitiya, N. Anousheh, S. Marru, G. C. Fox, and V. Jadhao, “Ions in nanoconfinement,” 2017, nanoHUB. [Online]. Available: https://nanohub.org/resources/nanoconfinement
 [57] G. C. Cawley and N. L. Talbot, “On overfitting in model selection and subsequent selection bias in performance evaluation,” Journal of Machine Learning Research, vol. 11, no. Jul, pp. 2079–2107, 2010.
 [58] M. Hunt, B. Haley, M. McLennan, M. Koslowski, J. Murthy, and A. Strachan, “Puq: A code for nonintrusive uncertainty propagation in computer simulations,” Computer Physics Communications, vol. 194, pp. 97–107, 2015.
 [59] J. Nossent, P. Elsen, and W. Bauwens, “Sobol’sensitivity analysis of a complex environmental model,” Environmental Modelling & Software, vol. 26, no. 12, pp. 1515–1525, 2011.
 [60] S. IM, “Sensitivity estimates for nonlinear mathematical models,” Math. Model. Comput. Exp, vol. 1, no. 4, pp. 407–414, 1993.
 [61] T. Homma and A. Saltelli, “Importance measures in global sensitivity analysis of nonlinear models,” Reliability Engineering & System Safety, vol. 52, no. 1, pp. 1–17, 1996.
Comments
There are no comments yet.