Atomic radius and charge parameter uncertainty in biomolecular solvation energy calculations

TL;DR

This paper introduces a novel method using surrogate models and generalized polynomial chaos to quantify uncertainty in atomic radii and charges, improving the reliability of solvation energy calculations in biomolecular modeling.

Contribution

It presents a new approach combining surrogate models, least-squares fitting, and compressed sensing to efficiently quantify parameter uncertainty in implicit solvation energy calculations.

Findings

Surrogate models effectively quantify uncertainty in small molecule solvation energies.

The method enhances understanding of parameter sensitivity in force field development.

Application demonstrates improved uncertainty estimation in biomolecular simulations.

Abstract

Atomic radii and charges are two major parameters used in implicit solvent electrostatics and energy calculations. The optimization problem for charges and radii is under-determined, leading to uncertainty in the values of these parameters and in the results of solvation energy calculations using these parameters. This paper presents a new method for quantifying this uncertainty in implicit solvation calculations of small molecules using surrogate models based on generalized polynomial chaos (gPC) expansions. There are relatively few atom types used to specify radii parameters in implicit solvation calculations; therefore, surrogate models for these low-dimensional spaces could be constructed using least-squares fitting. However, there are many more types of atomic charges; therefore, construction of surrogate models for the charge parameter space requires compressed sensing combined with an iterative rotation method to enhance problem sparsity. We demonstrate the application of the method by presenting results for the uncertainties in small molecule solvation energies based on these approaches. The method presented in this paper is a promising approach for efficiently quantifying uncertainty in a wide range of force field parameterization problems, including those beyond continuum solvation calculations.The intent of this study is to provide a way for developers of implicit solvent model parameter sets to understand the sensitivity of their target properties (solvation energy) on underlying choices for solute radius and charge parameters.

Tables

Table 1. Table 1: List of 17 compounds from the SAMPL computational challenge for solvation energy prediction with solvation energy [ 16 ] and solvent accessible volume from APBS [ 37 ] .

Ind.	Compound	Solvation energy (kJ/mol)	Molecular volume (Å3)
1	glycerol triacetate	-36.99	215.80
2	benzyl bromide	-9.96	126.98
3	benzyl chloride	-8.08	124.66
4	$m$-bis(trifluoromethyl)benzene	4.48	290.91
5	$N,N$-dimethyl-$p$-methoxybenzamide	-46.07	188.44
6	$N,N-4$-trimethylbenzamide	-40.84	179.62
7	bis-$2$-chloroethyl ether	-17.70	121.45
8	$1,1$-diacetoxyethane	-20.79	148.50
9	1,1-diethoxyethane	-13.72	132.75
10	$1,4$-dioxane	-21.13	87.86
11	diethyl propanedioate	-25.10	165.44
12	dimethoxymethane	-12.26	81.90
13	ethylene glycol diacetate	-26.53	148.10
14	$1,2$-diethoxyethane	-14.81	132.91
15	diethyl sulfide	-6.49	108.15
16	phenyl formate	-15.98	126.25
17	imidazole	-41.05	67.36

Equations

$$-\nabla \cdot ϵ(\mathbf{x})\nabla \phi (\mathbf{x})$$

$$\phi (\mathbf{s})$$

$$\Delta G={\int }_{\Omega }\rho (\mathbf{x})\left(\phi (\mathbf{x})-{\phi }_0(\mathbf{x})\right)d\mathbf{x},$$

$$\Delta G=\sum _i^{N_A}{q}_i\left(\phi ({\mathbf{x}}_i)-{\phi }_0({\mathbf{x}}_i)\right).$$

$$\text{Cov}({\mathbf{x}}_i,{\mathbf{x}}_j)={\eta }_i{\eta }_j\exp \left(-\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j{\Vert }_2^p}{\theta }\right),$$

$$\mathbf{q}=⟨\mathbf{q}⟩+{\mathbf{L}}_c\mathbf{\gamma },$$

$$\mathbf{\sigma }=⟨\mathbf{\sigma }⟩+{\mathbf{L}}_r\mathbf{\zeta },$$

$$RMS{E}_1=\sqrt{\frac{{\sum }_{q=1}^{10000}{\left(\tilde{E}({\mathbf{\xi }}^q)-E^q\right)}^2}{{\sum }_{q=1}^{10000}(E^q{)}^2}},RMS{E}_2=\sqrt{\frac{{\sum }_{q=1}^{10000}{\left(\tilde{E}({\mathbf{\xi }}^q)-E^q\right)}^2}{10000}}.$$

Full text

1 Introduction

Implicit solvent models and their applications have been the subject of numerous previous reviews [1, 2, 3]. Such solvation models require the coordinates of the solute atoms as well as atomic charge distributions and a representation of the solute-solvent interface. Charges and interfaces are generally modeled through parameterized empirical representations; however, these parameterizations are often under-determined, leading to uncertainty in the resulting parameter sets [4, 5, 6]. The Poisson equation is a popular model for implicit solvent electrostatics and serves as a good example for exploring the influence of this uncertainty on properties such as molecular solvation energy [1, 2, 3, 7]. In this paper, we use the term “solvation energy” to refer to the energy returned from the Poisson equation and to emphasize that we are not sampling over solute conformational states for a true free energy. This is a partial differential equation for the electrostatic potential $\varphi:\Omega\mapsto\mathbb{R}$

[TABLE]

where $\Omega\subset\mathbb{R}^{3}$ is the problem domain, $\partial\Omega$ is the domain boundary, $\epsilon:\Omega\mapsto[1,\infty)$ is a dielectric coefficient, $\rho:\Omega\mapsto\mathbb{R}$ is the charge distribution, and $\varphi_{D}$ is a reference potential function (e.g., Coulomb’s law) used for the Dirichlet boundary condition. The dielectric coefficient $\epsilon$ is usually defined implicitly [8, 9, 10, 11] with respect to the solute atomic radii $\{\sigma_{i}\}$ and solvent properties such that the coefficient reaches two limiting constant values: $\epsilon_{u}$ inside the solute and $\epsilon_{v}$ away from the solute in bulk solvent. The solvation energy is calculated by

[TABLE]

where $\varphi$ is the Poisson equation solution for the system with a bulk value of $\epsilon$ corresponding to the solvent of interest and $\varphi_{0}$ is the solution for the system with a bulk value of $\epsilon$ corresponding to a vacuum. For atomic monopoles, the solute charge distribution has the (numerically unfortunate) form $\rho(\bm{x})=\sum_{i}^{N_{A}}q_{i}\delta(\bm{x}-\bm{x}_{i})$ for $N_{A}$ solute atoms with positions $\{\bm{x}_{i}\}$ and charges $q_{i}$. The $\delta$ terms are formally defined as Dirac delta functionals but usually approximated by functions with finite support (e.g., when projected onto a grid or finite element basis). The delta functional approximation leads to a simplified form for the solvation energy in Eq. 3,

[TABLE]

Charges and interfaces in implicit solvent representations are generally modeled through parameterized empirical representations; however, these parameterizations are often under-determined leading to uncertainty in the resulting parameter sets [4, 5, 6]. For example, atomic charge models are designed to approximate the “true” vacuum electrostatic potential due to quantum mechanical electron and nuclei charge distributions. While quantum mechanical charge distributions can be incorporated directly in implicit solvent models [12, 13], atomic point charge distributions are generally used [2]. These point charges can include inducible and fixed multipoles [14, 15] but monopoles are the most common form. For the purposes of assigning charges, atoms are grouped into sets based on molecular connectivity and environment [16]. The charge values for atoms in these sets are usually determined by numerical fitting to quantum mechanical vacuum electrostatic potentials. Such charge optimization is ill-posed and fitting requires careful choice of the objective function and regularization constraints [17, 18, 19, 20]. While sophisticated fitting procedures have been developed, significant information reduction occurs in the transformation of the continuous quantum mechanical electron density into a discrete set of atomic point charges.

Solute-solvent interface models are much more empirical than the charge distribution models; the definition of a solvent “interface” is imprecise at length scales comparable to the size of water molecules. Therefore, such models are generally developed to represent a reasonable description of the solute geometry while also optimizing agreement with experimental quantities such as solvation energy. A large number of solute-solvent interface models exist, including van der Waals [11], solvent-accessible [21], solvent-excluded (or Connolly) [22], Gaussian-based [23], spline-based [24], and differential geometry surfaces [10, 25, 26, 27, 28, 29]. All of these interface models represent atoms as spheres and require information about the radii of these spheres. These radii are generally assigned to sets of atoms based on their “type” as determined by the local molecular connectivity. Unlike atomic charges, there are relatively few sets of atom types used to assign radii [16, 30]. These radii parameters are determined by optimization of properties such as solvation energy against experimental data [16, 30]. Additionally, many of these models also require information about solvent characteristics, generally in the form of a solvent radius, characteristic solvent length scales, or bulk solvent pressure/surface tension properties.

The intent of this study is to provide a way for developers of implicit solvent model parameter sets to understand the sensitivity of their target properties (solvation energy) on underlying choices for solute radius and charge parameters. In the present work, we present a new method to quantify the uncertainty in solvation energy calculated by the Poisson equation and induced by the uncertainty of the input radii and charge parameters. In particular, we construct two surrogate (or statistical regression) models of the solvation energy in terms of the radii and the atomic charges, respectively. These surrogate models enable us to estimate the solvation energy with different input parameters quickly and to evaluate the statistical information of the target properties (e.g., probability density function) efficiently. We model the input parameters as independent (i.i.d.) Gaussian random variables with different means and standard deviations; however, other probability distributions can also be used. To construct the surrogate of the Poisson model, we use a generalized polynomial chaos (gPC) [31, 32] expansion to represent the dependence of the solvation energy on uncertain parameters such as the atomic charge and radii. The efficacy of the gPC method for elliptic problems such as the Poisson equation has been extensively studied with robust results for its efficiency and accuracy [33, 34]. This approach is straightforward to apply to the relatively low-dimensional parameter sets. However, the main challenge of applying this method to implicit solvent calculation parameter uncertainty is the high-dimensionality of parameter sets (especially the atomic charges): the surrogate models require more basis functions and, therefore, more expansion coefficients need to be identified. To address this challenge, we adopt a compressive sensing method combined with the rotation-based sparsity-enhancing method first proposed by Lei et al. [35] and extended by Yang et al. [36], which enable us to construct the surrogate with relatively few sample outputs of the numerical Poisson solver.

2 Methods

We demonstrated the framework using a test set of 17 compounds from the SAMPL computational challenge for solvation energy prediction [16] (see Table 1). This set was chosen to demonstrate the uncertainty quantification framework on several different molecules; however, it was not chosen to calculate statistics over this small set.

We use this subset of the SAMPL data to demonstrate the use of our method to quantify uncertainty in solvation energy due to implicit solvent parameter uncertainty.

2.1 Uncertain parameters

Many parameterization approaches for atomic charge use ESP (electrostatic potential) [38] or related methods (e.g., RESP [19]). These methods optimize atomic charges by least-squares fitting of the charges’ Coulombic potential to the electrostatic potential obtained from quantum mechanical calculations. This under-determined optimization is performed subject to various constraints, including the requirement that the atomic charges sum to the integer formal charge of the molecule. More specifically, the calculated ESP $\hat{V}_{i}$ at the $i$-th grid point is the electrostatic potential given by Coulomb’s law summed over the charge $q_{j}$ at the centers of the $j$-th atoms. Least-squares fitting is performed by minimizing $\sum_{i}(V_{i}-\hat{V}_{i})^{2}$ with constraints, where $V_{i}$ is the electrostatic potential computed by ab initio calculations. Least-squares fitting implies a Gaussian noise model wherein the atomic charges $q_{j}$ can be modeled as Gaussian random variables as done in this study.

In the present work, we modeled the uncertainty in atomic charges by considering atomic charges obtained by 11 different approaches: AM1BCC [39], CHELP [40], CHELPG [41], CM2 [42], ESPMK [38], Gasteiger [43], PCMESP [44], Qeq [45], RESP [19], MMFF94 [46], Mulliken [47]. The Hartree-Fock method and the 6-31G*basis set were used to optimize molecular geometries. The methods we selected here are popular for implicit solvation models and all-atom approaches. Although many of these charge methods are used in all-atom simulations, implicit solvent models have been used with several of them, including RESP and ESP(MK) [48, 30], AM1-BCC [48], Mulliken [49], CHELPG [50, 49], Gasteiger [51], Qeq [52], etc.

We have assumed that the variation of atomic charges across different methods can be modeled by a Gaussian random field with covariance kernel

[TABLE]

where $\eta_{i}$ is the standard deviation of the $i$-th atomic charge, $\bm{x}_{i}$ is the position of the $i$-th atom, and $0<p<2$. The least-squares nature of most charge fitting methods makes Gaussian variables a natural choice; however, other probability distributions can also be used. We used atomic charges from 11 different methods to estimate $\eta_{i}$ and then used the maximum likelihood estimate (MLE) method to estimate $\theta$ and $p$. Since the sum of $N_{A}$ charges in a molecule is constrained (to its formal molecular charge $Q\in\mathbb{Z}$), we modeled the Gaussian random field with $N_{A}-1$ atoms by removing the last hydrogen in the PDB file. Additionally, we use symmetry in the molecular structure to reduce the number of independent atomic charge types before applying the MLE to identify the random field. For example, in a benzene, there is only one type of carbon and one type of hydrogen due to the symmetry of this molecule. Therefore, we considered the charges of its atoms as a Gaussian random field with only two entries instead of 12 ones (the total number of atoms in benzene).

After obtaining the covariance matrix by integrating across methods, we represented the atomic charge as

[TABLE]

where $\bm{q}=(q_{1},q_{2},\cdots,q_{{}_{N_{A}-1}})$ are the atomic charges, $\left\langle{\bm{q}}\right\rangle$ is the mean of $\bm{q}$ estimated from the 11 different charge values, $\bm{\gamma}=(\gamma_{1},\gamma_{2},\cdots,\gamma_{{}_{N_{A}-1}})$ are i.i.d. zero-mean unit-variance Gaussian random variables, and $\bm{L}_{c}$ is a lower triangular matrix from the Cholesky decomposition of the covariance matrix (Eq. 5). We note that for the atoms in the test set used in the present work, the covariance matrices of these random field are almost diagonal: the off-diagonal entries are smaller than $10^{-12}$. This suggests the correlation between atomic charges is effectively removed during their symmetry-based grouping. The atomic charge for the remaining atom is obtained by summation of the other random charge variables based on the constraint $q_{i}=Q-\sum_{j\neq i}^{N_{A}}q_{j}$.

Similarly, we used multiple force fields (ZAP-9 [16], OPLSAA [53], Bondi [54] and PARSE [30]) to model uncertainty in the radii parameters in the same manner. Although radii are non-negative, we did not explicitly impose constraints on the radii. After obtaining the covariance matrix, we represented the radii as

[TABLE]

where $\bm{\sigma}=(\sigma_{1},\cdots,\sigma_{N_{A}})$, $\sigma_{i}$ is the radius of atom (type) $i$, $\bm{\zeta}=(\zeta_{1},\cdots,\zeta_{N_{A}})$ are independent zero-mean unit-variance Gaussian random variable and $\bm{L}_{r}$ is a lower triangular matrix from the Cholesky decomposition of the covariance matrix. The small number of radii sets makes the selection of a probability distribution somewhat arbitrary. We have assumed that the radii follow a Gaussian distribution; however, other probability distributions can also be used. We note that the standard deviations here are smaller than $10\%$ of the mean values which implies very low probabilities for unphysical negative radii values. Therefore, by employing truncated Gaussian random variables within 4 standard deviations (capturing more than 99.99% of the probability), we guaranteed that the radii are always positive and that the distributions of the truncated Gaussian variables were almost identical to the original Gaussian variates. We note that with this setting, no model assigns zero radius to protons or other atoms.

Although we use $\gamma$ and $\zeta$ to denote the random variables used for modeling the uncertainties in $q_{j}$ and $\sigma_{j}$, in what follows, we still use $\bm{\xi}=(\xi_{1},\xi_{2},\cdots)$ to denote general uncertain inputs when introducing the algorithm and reporting results.

2.2 Solvation energy surrogate models

We used generalized polynomial chaos (gPC) expansions as surrogate models for the solvation energy. The goal of surrogate construction is to estimate the variations in quantities of interest, such as solvation energy, much more efficiently than solving the original problem, such as solving the Poisson equation. The details for these expansions are provided in Supporting Material.

2.3 Poisson equation solver

We used the Adaptive Poisson-Boltzmann Solver (APBS) [37] to solve the Poisson equation for solvation energies. Poisson calculations were performed with the finite difference solver using $97^{3}$ grids focused from a 25 Å to a 13 Å cubic domain. Charges were discretized onto the grids using linear interpolation. Boundary conditions were assigned using a sum of Coloumb potentials. The molecular interior and solvent were assigned dielectric values of 2.0 and 78.0, respectively. The solute-solvent boundary was defined using a “Connolly” molecular surface [55]. Energies were calculated using the standard approach for Poisson-Boltzmann calculations [56, 57].

3 Results and discussion

For each test case, we used Monte Carlo simulations to sample the parameter probability distributions and generate $10,000$ samples of the input parameters $\bm{\xi}^{q}$ and then solved PB equation using APBS to obtain output samples of the solvation energy $E^{q}=E(\bm{\xi}^{q})$. We used these outputs as ground-truth reference solutions to examine the performance of the surrogate models; these outputs will be referred to as “reference” in the remainder of this paper. More precisely, given a surrogate model $\tilde{E}$, we use two different root-mean-squared error (RMSE) measures to examine its accuracy:

[TABLE]

We also use box-whisker plots to demonstrate the statistics. The line in the middle is the median of 16 molecules, the tops and bottoms of the boxes are 25th and 75th percentiles, and the whisker plots cover more than $99\%$ probability.

3.1 Influence of radii uncertainties on solvation energies

We investigated the effect of the uncertainties in the radii with fixed atomic charges obtained from AM1-BCC [39]. As an example, there are eight different sets of radii for $N,N$-dimethyl-$p$-methoxybenzamide across the ZAP-9, Bondi, OPLSAA, and PARSE parameter sets, as shown in the support material. We modeled the solvation as a function of eight i.i.d. Gaussian random variables. We constructed gPC surrogate models with multi-variate normalized Hermite polynomials up to third order. The surrogate model consisted of $C_{8+4}^{4}=495$ basis functions. Figure 1 (a) presents the RMSE obtained by our method with respect to different numbers of samples $E^{q}$. Figure 1 (b) compares the solvation energy probability distribution function (PDF) obtained by our method and the reference solutions. The numerical results are obtained by constructing the surrogate model with the $36$ output samples first, then sampling the surrogate model $10,000$ times with random samples to estimate the PDF. The reference solution is computed from the $10,000$ outputs of $E^{q}$.

We performed the same analysis for all the molecules in the test set and present the results in Figure 2. For most molecules, we can build an accurate surrogate model (RMSE$<0.05$) for the solvation energy with only a few samples (less than $40$) of the input parameters. However, $m$-bis-trifluoromethylbenzene (TFMB) required significantly more samples. In particular, the RMSE for the TFMB solvation energy surrogate model was close to $0.15$ with $40$ samples and required $100$ samples to reduce the RMSE to less than $5\%$. This variability arises from the radius of fluorine: in the ZAP force field it is 2.4 Å; however, it is only $\sim 1.4$ Å for the other force fields. Hence, the standard deviation of this radius is around $25\%$ of the mean and fluorine requires more terms in the surrogate model for an accurate description and therefore more samples to parameterize those terms.

The influences of the uncertainties in the input radii on the solvation energy for each molecule are demonstrated in box-whisker plots in Figure 3. The experiment results are presented for comparison. We note that some experiments results are “outliers” of the box-whisker plots, this is because that the atomic charges are computed from AM1BCC for the purpose of fixing the atomic charges and it does not guarantee that the computed solvation energy is sufficiently close to the experiment results. For example, for the $m$-bis(trifluoromethyl)benzene AM1BCC charges yield negative solvation energy while the experiment result is positive.

3.2 Influence of atomic charge uncertainties on solvation energies

We also examined the influence of charge perturbation for solvation energy calculations with fixed radii (ZAP-9). As an example, there are 14 different types of atoms in $N,N$-dimethyl-$p$-methoxybenzamide as shown in Supporting Material. We note that we model the surrogate with 13 inputs due to the constraint on the summation of the charges. The mean and standard deviation are computed from the results of 11 different charge fitting approaches. We used no more than 3000 multi-variate normalized Hermite polynomials (up to fourth order) in the gPC surrogate model for $E_{g}$ for all the molecules. We use $N,N$-dimethyl-$p$-methoxybenzamide as an example. Figure 4 (a) presents the RMSE obtained by our method with respect to different numbers of samples $E^{q}$. It illustrates that 300 output samples are needed to reduce the RMSE to less than $5\%$. Figure 4 (b) compares the PDF obtained by our method and the reference solution. The numerical results are obtained by constructing the surrogate model with the $300$ output samples first, then sampling the surrogate model $10,000$ times with random samples to estimate the PDF. The reference solution is computed from the $10,000$ outputs of $E^{q}$.

The influences of the uncertainties in the input atomic charges on the solvation energy for each molecule are demonstrated in Figure 5. For most molecules, the experiment results lie in the whisker plots and some of them are in the box.

We also present the number of output samples needed to construct a surrogate with RMSE less than $5\%$ with respect to the number of atom types in Figure 6.

3.3 Combined influence of radius and atomic charge uncertainties

We chose the charge methods and radii based on their popularity in the implicit solvent community. Not all of the radii and charges examined in this study would be expected to give accurate answers when used together. We could have chosen a more constrained set; however, we chose this diverse set as a more challenging example to test our method to illustrate our approach across significant variation in parameter values. Comparing the PDFs in Figures 1 (b) and 4 (b), we notice that the uncertainty in the solvation energy induced by the atomic charges is stronger than that induced by the radii. A similar observation has been made previously by Chakavorty et al. [58] who also noted that conformation introduces another important source of uncertainty across force fields. We have investigated the influence of conformational uncertainty on solvation in a previous paper using similar methods [35]; however, combining parameter and conformational uncertainty is outside the scope of this manuscript.

The atomic charges vary significantly across different methods while the variation in the radii is much smaller. To understand the combined influence of charges and radii on solvation energies, we modeled the correlated uncertainties for these two types of parameters can be modeled with i.i.d. Gaussian random variables. We use $N,N$-dimethyl-$p$-methoxybenzamide as an example. 480 output samples are needed to reduce the RMSE to less than $5\%$. Figure 7 (a) presents the RMSE obtained by our method with respect to different numbers of samples $E^{q}$. Figure 7 (b) compares the PDF obtained by our method and the reference solution. The numerical results are obtained by constructing the surrogate model from the $480$ output samples and then sampling the surrogate model $10,000$ times with random samples to estimate the PDF. The reference solution is computed from the $10,000$ outputs of $E^{q}$.

Not surprisingly, the number of output samples needed to construct an accurate surrogate increases as we take into account both uncertainties in the charges and radii. The shape of the solvation energy changes PDF also slightly as the radii variation of the radii across different methods are much smaller than charge variations.

The influences of the uncertainties in the input atomic charges on the solvation energy for each molecule are demonstrated in Figure 8. This figure is similar to Figure 5 since the uncertainties in the atomic charges dominate the results.

Figure 9 shows the number of output samples needed to construct a surrogate with less than $5\%$ RMSE for all 17 molecules in the test set.

This figure illustrates the approximately quadratic scaling with the respect to the number of atom types in the molecule.

4 Conclusions

We have developed a new method for quantifying the uncertainty associated with parameterization of implicit solvent models. In particular, we used a newly developed extension of compressive sensing method to construct surrogate models of solvation energy based on gPC expansions. These surrogate models allow us to efficiently and accurately estimate the variation in solvation energy due to uncertainty in charge and radius parameters. In this initial work, we used statistical distributions for radius and charge variation based on the observed differences in the parameter sets. However, in future studies, it may be useful to use the uncertainty quantification approach presented here with more physically motivated models that address the underlying uncertainties in determining charge and radius parameters. Our results demonstrate that for the data sets used in the present work, the variation of radii across different approaches are small. On the other hand, the variations of the atomic charges obtained by different methods are much larger, so that the number of output samples needed for accurate UQ analysis requires are much larger, growing quadratically with respect to the number of atom types. This framework can be applied to estimate the statistics (e.g., mean, variance), PDF, confidence interval, Chernoff-like bounds [59], etc. of solvation computing and other chemical computing when the inputs are uncertain. The current study focused on uncertainty in solute charges and radii; however, this framework could also be applied to other solvation model characteristics such as dielectric coefficient, solvent radius, and biomolecular surface definition. Likewise, this approach could also be used for quantities of interest other than solvation energy; e.g., dipole moments, titration states, etc.

In the future, we anticipate that this approach could be used for a much wider range of force field parameterization activities, including both coarse-grained and atomistic representations of biomolecules. Uncertainty quantification methods have begun to be used in force field parameterization of simple alkane systems [60]; this paper demonstrates the ability to extend the methods to higher-dimensional systems with more diversity of atom types. Application of these methods offer the benefit of efficiently characterizing parameter space and understanding the impact of parameter variation on quantities of interest. Additionally, the iterative method we used in the present work is very suitable for this type of problem, as the accuracy of the surrogate models are improved significant after iterations. Especially, the error of the surrogate models for the atomic charge induced uncertainties are reduced by $40\%\sim 50\%$ compared with the standard compressive sensing method. Also, there is significant room for development in the numerical methods. For example, the sparsity-enhancing approaches can be combined with other techniques including improved sampling strategies [61, 62], adaptive basis selection [63, 64], and advanced optimization methods [65, 66]. These approaches improve the accuracy of the compressive sensing method from different aspects. As such, they will help to reduce the number of expensive simulations or quantum mechanics calculations needed for constructing accurate surrogates.

Acknowledgments

This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific-Computing Research as part of the Collaboratory on Mathematics for Mesoscopic Modeling of Materials (CM4) and by NIH grant GM069702. Pacific Northwest National Laboratory is operated by Battelle for the DOE under Contract DE-AC05-76RL01830.