Adaptive selection and validation of models of complex systems in the presence of uncertainty

This paper describes versions of OPAL, the Occam-Plausibility Algorithm (Farrell et al. in J Comput Phys 295:189–208, 2015) in which the use of Bayesian model plausibilities is replaced with information-theoretic methods, such as the Akaike information criterion and the Bayesian information criterion. Applications to complex systems of coarse-grained molecular models approximating atomistic models of polyethylene materials are described. All of these model selection methods take into account uncertainties in the model, the observational data, the model parameters, and the predicted quantities of interest. A comparison of the models chosen by Bayesian model selection criteria and those chosen by the information-theoretic criteria is given.

ciple of parsimony of explanations"). So-called Occam factors as a measure of the relative value of one model over another are discussed in, for example, Jaynes [15], Loredo [23], Wolpert [27,28], and elsewhere as a form of a Bayes factor weighted by maximum likelihoods of the respective models. Our own version, embedded in OPAL, provides for not only a parsimonious approach, but, importantly, an approach for addressing model inadequacy and determining model validity relative to estimates of the accuracy with which the model predicts representations of quantities of interest. OPAL is based on Bayesian model plausibilities, but also involves partitioning models into "Occam Categories," which are derived from a measure of simplicity based on the number of parameters in a model.
In the art and science of model validation, much depends upon how one determines that a model "adequately represents" observational data. Such a determination requires a notion of adequacy, i.e., a measure of accuracy with which a model can predict specific data, and a tolerance that must be met in order to deem a model sufficiently accurate. A famous quote, also attributed to Box [6], is "all models are wrong, but some are useful." Validation processes aim to judge which models are useful in predicting specific events in physical systems.
We remark that many studies have been performed on methods of model selection in statistics and biological literature, good examples being the work of Posada and Buckley [25] on model selection and averaging in phylogenetics; the work of Gelman, Hwang, and Vehtari comparing Akaike information criterion (AIC), deviance information criterion (DIC), and Watanabe-Akaike information criterion (WAIC) approximations of crossvalidation [12]; the book of Burnham and Anderson [8] on ecological models; and the book of Konishi and Kitagawa [22] on information-theoretic methods in statistical modeling. These studies do not address the fundamental issue that the best model in a set of models, no matter how "goodness" of the model is measured, may be completely unacceptable for the predictive purpose at hand; i.e., the best model may be invalid. The approaches described in the present work address both relative model quality and validity.
In the present paper, we examine and compare alternative forms of OPAL in which different methods of model selection are employed. In particular, we depart from a fully Bayesian approach and explore the frequentist, information-theoretic methods embodied in the Akaike information criterion (AIC). Introduced in the 1970s [2][3][4], AIC has been studied and used predominantly in areas of ecological and biological sciences, as discussed, for example, in [8]. Variations and extensions have also been introduced to confront computational complications that may arise in cases of limited observational data; see, e.g., [14,22,26]. We give further details on AIC, Bayesian plausibilities, and other methods and compare validation results and predictions using each approach. We continue to focus on the difficult problem of selection and validation of coarse-grained models of atomistic systems, as it exhibits all of the challenges of model validation and quantification of uncertainties in predictions.
Following this Introduction, we review a number of basic concepts that are fundamental to predictive science and, particularly, model validation. We review general methods of model selection in Sect. 3 and describe OPAL in Sect. 4. Applications to coarse graining of models is taken up in Sect. 5 and conclusions are collected in a final section.

Preliminaries
To establish the setting for this exposition, we review briefly some basic concepts and notations relevant to predictive modeling following [24]. It is understood in the present work that the term model refers to a mathematical model, a collection of mathematical constructions put forth as mathematical abstractions of systems, particularly physical or engineered systems, but social and economic systems could conceivably be considered as well. It is useful to regard a model as characterizing an abstract mathematical problem, such as: Given θ ∈ ⊂ R k and S, find u(θ, S) ∈ V such that where A is a collection of mathematical operators, θ is a vector of model parameters taken from a parameter space (assumed finite dimensional here), S is the scenario in which the model is implemented, and u(θ, S) is the solution for given θ and S belonging to a space V of trial functions. For any S, (1) is referred to as the forward problem, as the model extrapolates information forward into a solution u(θ, S) from which predictions of features of the system (or specific events) are derived. The concept of a scenario is important in the science and technology of model validation. Mathematically, a scenario is viewed as a set of parameter-independent features of the model that can generally be specified exactly, such as the domain of the solution u(θ, S) or certain boundary and initial data, the idea being that the same model can be used in several different scenarios. The term scenario is used to refer to both the actual physical environment, in which experimental data are collected, and the computational environment, in which the reality to be predicted by the model resides.
In theory, the mathematical model is selected to solve the forward problem in the full prediction scenario S p , and the solution u(θ, S) is then used to compute specific quantities of interest (QoIs). The QoIs are specific realities targeted in the prediction process. Mathematically, they are usually characterized as functionals Q on the space V of functions in which the solution u(θ, S p ) resides: Physically (and philosophically), the QoIs are not observables, as model predictions may be of events in systems that do not (yet) exist, such as engineering designs. In statistics, they are sometimes identified as "out-of-sample" data (e.g. [12]), extrapolations outside realm of measurable observations. In processes of calibration and validation of models, other scenarios are considered. At a primitive level, calibration scenarios S c are considered that involve unit tests on components of a model. They are designed to update prior information on model parameters by matching model predictions with experimental calibration data y c . Multiple calibration scenarios may be considered, each involving different model parameters that could be subsets of those appearing in the prediction scenario. The fundamental issue of model validation is the process of assessing the validity of the model in question as a means to predict the QoI with acceptable accuracy. This involves the design of validation experiments on subsystem models in validation scenarios S v , designed to compare model predictions with validation observable data y v . The challenges of validation are to design experiments that deliver observational data adequately representing the QoI, to test the validity of hypotheses made in developing the model that may not be fully trusted, to choose an appropriate metric to measure the accuracy with which the model predicts QoI-informed data y c , and to select a tolerance γ tol , of the error that the modeler is willing to accept in order to declare the model "valid" (or not invalid). Once a model is determined to be valid (through this subjective process), the model parameters of the valid model (or their statistical representation by appropriate probability density functions) are introduced into the model implemented in the full prediction scenario, the forward problem is solved, and the QoI is evaluated.
Several remarks should be made at this point. Firstly, as noted earlier, and famously noted by Box [6], all models of physical reality are imperfect. The goal of predictive computational science is to determine whether model predictions are "close enough" to reality to use in predictions of events to contribute to scientific knowledge or as a basis for important decisions. Next, the validity of a model depends upon the QoI to be predicted; a model "valid" for one QoI may be invalid for another. It is emphasized that the notion of a valid model can be highly subjective, involving the design of an experiment to mimic a non-observable QoI, the choice of a metric, and the choice of a tolerance for acceptability. Furthermore, predictions are made in the presence of many uncertainties in model parameters, in data (y c , y v ), and in selecting the model itself. Validation methods may or may not address these uncertainties. The QoI is generally a random number or variable. An important challenge is to quantify the uncertainty in the prediction. In addition, a computational model is generally derived from the mathematical model to render it to a form that can be processed by a computer. The discretization of the model, of course, introduces additional errors in the prediction. While not addressed in the current study, this subject is taken up in earlier work [1]. Finally, not all of the parameters of a model necessarily influence the QoI for a particular prediction scenario S p . Effective methods of measuring the sensitivity of predictions to choices of model parameters and estimates of parameter sensitivity can lead to elimination of models that do not significantly inform the QoI, resulting in a substantial reduction in the complexity of the model selection process.

Model selection
The development of a rigorous basis for selecting the "best" model among a set of possible models has been a goal of some modelers, particularly in statistics, for decades. Various measures to assess the quality of one model relative to another have foundations in Bayesian arguments, such as Bayes' factors and Occam factors, as well as information-theoretic arguments derived from frequentist statistics and maximum likelihood approaches. All of the methods of interest involve calculations designed to assess how well model predictions using parametric models agree with observational data or how closely the parametric model can approximate a probability distribution representing the "truth," i.e., the true reality.
In the Bayesian setting, the idea of model evidence and posterior model plausibilities is extremely powerful. The notion of posterior plausibilities is mentioned in the 1981 paper of Chow [9], who attributes the idea to Jeffreys' treatments of probability theory [16]. It was certainly known to Schwarz [26], who developed easily implemented approximations to model evidence that lead to the Bayesian information criterion (BIC) in analogy to information-theoretic approaches. More recently, the use of such Bayesian probability approaches for model selection were advocated by Beck and Yuen [5], Hawkins-Daarud et al. [13], and Farrell et al. [10,11]; see also [24].
On the information-theoretic side, the work of Akaike [2,4] leading to the Akaike information criteria or its various generalizations [14] are perhaps best known in the domain of frequentist statistics. An account of information-based model selection criteria, including extensions to "generalized information criteria" (GIC), is given in the book of Konishi and Kitagawa [22].
We explore the case in which we have many models from which to choose, and the goal is to identify the model with the most potential for predicting target quantities of interest. Consider the set M of m parametric model classes, each with its own parameter space i from which the parameter vector θ i may be chosen. Each model P i (θ i ) is presumed to be developed on the basis of physical and empirical laws described by mathematical constructions. Assume also that we have acquired a set of observational data y = {y 1 , y 2 , . . . , y n } ∈ Y S where Y S is a space of observational data accessible in scenario S.

Bayesian posterior plausibilities
Suppose that for each model, a prior probability density π(θ i |P i , M) is specified. We simply rewrite Bayes' rule, acknowledging that we have additional conditional information; namely, that the model P j (θ j ) is known to belong to the set M: The key term is the model evidence, the denominator on the right-hand side of (4). It is the marginalization of the numerator with respect to the parameters: The model evidence can be interpreted as a new likelihood function for a discrete version of Bayes' rule over the set M of models. Its posterior, for each model P j , denoted ρ j , is the posterior model plausibility, Choosing the denominator to normalize the set of discrete plausibilities, we have The model (or models) in M with the largest plausibility (or plausibilities) ρ k such that ρ k ≥ ρ j , 1 ≤ j ≤ m, is deemed the most plausible (the "best") model in the set M for given data y. The prior π(P j |M) may be set equal to 1/m if, initially, all m models are regarded as equally plausible. Otherwise, prior experience may be called upon to weigh one model over another.

The Akaike information criterion
As in the Bayesian setting, each model has its own likelihood distribution π j (y|θ j ). We drop the dependence on M and replace the dependence on P j with the subscript π j for the moment to simplify notation. Note that each likelihood captures the probability that the model P j is able to reproduce the observed data y.
It is customary and convenient in mathematical statistics, particularly in frequentist statistics, to assume the existence of the "truth" or full reality embodied in a probability density f . The expected value with respect to the truth f of the log-likelihood has the following important property: Denote by θ * the vector such that where D KL is the Kullback-Leibler divergence, also called the the relative entropy or the information loss between the truth f and the model, represented by the likelihood π(y|θ), That is, the parameter vector θ * that maximizes the log-likelihood with respect to f for given data y minimizes information lost in approximating the truth with the model. For a set M of m models, each with likelihood π j (y|θ j ), the averaged D KL over the truth leads to the following measure of information loss for each model: x and y being independent samples taken from the same distribution for scenario S and θ j (y) is the maximum likelihood estimate (MLE) for model P j , found using data y. As in (8), the expectations E x and E y are taken with respect to the truth f . Obviously, since the truth f is not known, the measures A j cannot be evaluated. However, under certain assumptions on smoothness of the likelihood function and asymptotic behavior of the likelihood as the number n of samples becomes larger, the quantities A j can be approximated by the Akaike information criterion, where K j is the number of parameters of model P j . The model with a lower AIC indicates a model with a lower D KL -distance to reality and is, therefore, a "better" model. Therefore, in general, we seek model P j ∈ M such that A derivation of (11) is given in [8].
An accepted alternative, the Bayesian information criterion (occasionally referred to as SIC, the Schwarz information criterion), is given by [26], Also, a "second-order" AIC for small sample sizes has been proposed by Hurvich and Tsai [14] which replaced AIC j of [7] by n being the sample size. Burnham and Andersen [8] comment that "unless the sample size is large with respect to the number of estimated parameters, use of AIC cj is recommended" over AIC j .
The individual values of AIC j , BIC j , or AICc j are generally not interpretable; however, the relative values are informative. To this end, and recalling that the minimum AIC j (or BIC j or AICc j ) value indicates to the "best" model, we can calculate the AIC differences, The larger j is, the less plausible it is that P j is the best model for minimizing the information lost in using model P j over the "truth" model. Recalling that the evidence π(y|P j , M) of (5) can be regarded as the likelihood of the model P j given the data y, an information-theoretic version of this likelihood can be put forth that is related to the AIC differences of the form, Then the Akaike version of model plausibilities is, in analogy with (6), the so-called Akaike weights, The Akaike weights w i are akin to Bayesian model plausibilities and can be used as a model selection criterion, and the model with weight w i closest to unity is deemed the best.

OPAL
As previously stated, OPAL is an algorithm designed to systematically select the simplest valid model. In the version advocated in [11], the simplicity was defined by the number of parameters such that the simplest model has the fewest parameters among a set of models, and validity was established by passing a validation criterion. It is clear from (11), is identified, each with parameters θ i belonging to an appropriate parameter space A parameter sensitivity analysis is performed to assess the sensitivity of a model output function Y (θ) on perturbations in model parameters. Those models with parameters not appreciably affecting the output are eliminated, yielding a reduced setM of models, M = P 1 (θ 1 ),P 2 (θ 2 ), . . . ,P l (θ l ) , l ≤ m.
3. The models surviving Step 2 are partitioned into "Occam Categories" according to their complexity. Those with the fewest parameters, for example, are put in Category 1, those with the next highest number of parameters in Category 2, and so forth. 4. Models in the set M * of Category 1 are calibrated in calibration experiments involving calibration data y c , yielding a calibrated set of Category 1 models, 5. The posterior Bayesian plausibilities ρ i of all models in M * are computed. Recall that these plausibilities depend explicitly (see, e.g., (6)) and implicitly (via the calibration process (4)) on the calibration data y c . Only the most plausible models with ρ i ≥ ρ j , 1 ≤ j ≤ m are retained.
6. An experimental validation scenario is constructed yielding validation observational data y v , and the most plausible model in Category 1 is used to compute a prediction of the observables y v ; if the difference between the observables and the prediction, measured in an appropriate metric or pseudo-metric, is within a preset tolerance γ tol , the model is deemed "valid." If not, one returns to Step 3 and repeats the process for the next category of models until a valid model is found. If no models of any category are deemed valid, one returns to Step 1 and enlarges the set M of possible model classes and then proceeds with the steps listed above. 7. Upon identifying a valid model, the forward problem is solved in the prediction scenario and the original QoI is computed, completing the prediction process.
All of these steps are designed to cope with uncertainties in the parameters, the observational data, and the target QoIs, all generally characterized by probability densities. The output Y (θ), when feasible, may be taken to be the QoI available in the full prediction scenario of the model. In regard to eliminating parameters with small sensitivity indices with respect to the output function Y (θ), it should be noted that the elimination of a parameter should be done only if 0 is in the domain of the parameter itself. Otherwise, another nominal value of the parameter must be chosen for the reduced model(s). It should be understood that among the set M of model classes, there may be many better models than that selected by OPAL, i.e. models in a higher Occam category, which produce predictions closer to the validation observations by some appropriate metric. OPAL is designed to uncover the simplest model (as measured by the number of parameters, for example) that satisfies a predefined validation criterion.
Step 4 in the OPAL algorithm is often the most computationally intensive, and it may be meaningful to consider other simpler methods of model selection when feasible. One goal of this study is to explore, through numerical experiments, the results of model validation when simpler methods, such as the AIC and BIC, are used instead of plausibilities for complex multi-parameter problems. Both the AIC and the BIC are derived using several simplifying approximations that involve truncation error and use of asymptotic estimates, and are not regarded to deliver model selections as accurate as plausibility measures.

Application to the selection of coarse-grained models of atomistic systems
One of the most complex challenges in model selection and validation occurs in the construction of coarse-grained (CG) models of atomistic systems-a standard approach in molecular dynamics simulations of chemical and biological systems. CG models are created by aggregating atoms together into representative groups. Interactions between these new groups are generally unknown and must be defined in terms of force potentials to characterize the mathematical representation of each CG model of the molecular system, the parameters of which should be determined following theories, ideas, and processes discussed earlier.
Traditionally, many interatomic potentials are represented by four types of interactions: bonded, angular, torsional (dihedral), and van der Waals. It is quite common to represent bonds and angles with a harmonic potential and the van der Waals interactions with a Lennard-Jones potential. In the OPLS functional form [17,18], torsional interactions are given a cosine expansion. The OPLS potential energy is therefore given by, In this equation, r is the configuration of the atoms (the vector of atomic coordinate vectors), the parameters k r,i and k θ ,i are spring stiffness parameters for the bonds and angles, respectively, r 0,i and θ 0,i are the equilibrium bond and angle values, and r i and θ i are the bond and angle distances at the given configuration r. The constants V n,i are the dihedral coefficients, and φ i are the torsional angles. For any pair of non-bonded atoms, the Lennard-Jones parameters are ij and σ ij , and r ij is the distance between them. An electrostatic interaction in the form of a Coulomb potential can also be added; however, in this example, the CG particles are assumed to be charge neutral.
In both the atomistic and coarse-grained systems, the potential energy drives the computational simulation. During these implementations, configurations are sampled and the corresponding potential energy is calculated. The probability density approximated using these samples will be used to the compute the validation metrics discussed later in this section.
Following [11], we consider as a representative example the problem of computing the potential energy of a polyethylene cube. To initialize OPAL, a set M of possible parametric models is identified. First, the AA-to-CG map is defined, shown in Fig. 1. Each CG particle, shown in red, is defined to represent two carbon atoms and their attached hydrogen atoms. In this example, the set M is created by tabulating the possible combinations of interactions shown in (21), as shown in Table 1. Further details regarding these types of interactions can be found in [10,11]. The interactions and parameters of each model are the same for those used in the Bayesian setting; therefore, the sensitivity analysis detailed and completed in [11] may be used here. As discussed in previous work, the potential energy in this particular application is insensitive to dihedral parameters. In some cases,  The Occam categories are determined by counting the number of parameters in the model insensitive parameters may be set to nominal, constant values. In the present example, the models are nested; thus, the models containing dihedral interactions are redundant and may therefore be eliminated from consideration. The remaining models are collected into the setM such thatM = {P 1 ,P 2 , . . . ,P 11 } = {P 1 , . . . , P 9 , P 11 , P 12 }. From Table 1, it can be seen that the lowest category contains those models which depend upon only three parameters. Thus, M * = {P * 1 , P * 2 , P * 3 , P * 4 } = {P 1 , P 2 , P 3 , P 4 }. The MLE for each of these models is determined via a quasi-Newton optimization scheme in which the starting point is the mean value of the parameters determined in an analysis of a simplified AA scenario. Specifically, these mean values are those used in the maximum entropy prior distributions in the Bayesian implementation. See appendix of [11] for complete details.
The calibration scenario consists of a single chain of polyethylene (C 80 H 162 ), simulated in a canonical ensemble. The calibration data y c is a vector of observed potential energy values. Once the MLE for each model in Category 1 is obtained, the AIC and Akaike weight for each model may be calculated according to (11), (15) and (17). In the present application, we find that For this set of models, it is clear that P * 1 = P 1 , which contains only bonded interactions, is considered the "best" model. This result agrees with the best model choice determined using Bayesian plausibilities for model selection [11].
Separate validation scenarios in which parameters are updated again are not typical in deterministic model development. However, for the purpose of following OPAL, we construct a validation test for the AIC-best model. We consider two chains of C 80 H 162 simulated in a canonical ensemble and the data y v is a set of potential energies. Using the calibration MLE as the starting point for the validation likelihood maximization, we update the MLE.
As validation metrics, we calculate the D KL between the AA and CG distributions produced in the validation scenario, as well as the normalized Euclidean distance between the AA and CG ensemble averages. That is, if u AA and u CG are the set of samples of the potential energy produced by the AA and CG models, respectively, in a given validation scenario, and We take as validation tolerances, for the normalized Euclidean distance and D KL metrics, respectively. Notationally, Q AA is the ensemble average of the observable, O(Q AA ) is its order of magnitude, and σ 2 AA is its variance.
For model P * 1 , the MLE parameters are updated using data from the validation scenario. These parameters are used in a forward implementation of polyethylene. Then, using (23) and (24), rendering P * 1 "valid" (not invalid). In a second validation scenario consisting of four chains of C 80 H 162 , the parameters, without another update, are used to produce a second u AA and u CG . Using the same validation tolerances γ 1,tol and γ 2,tol , both models are again deemed "not invalid" since Since this model has passed two validation tests, we have confidence to say that it may be used to predict the QoI in our prediction scenario. Although P * 1 has been vetted as a "valid" model, we shall move through another iteration of OPAL for the purposes of illustration. Tightening the validation criteria so that renders the Category 1 model P * 1 invalid. Following the OPAL algorithm, we move to the subset of models in Category 2, in which M * = {P * 1 , P * 2 , P * 3 , P * 4 , P * 5 } = {P 5 , P 6 , P 7 , P 8 , P 9 }. For each model, the MLE parameters and Akaike weights are calculated for each model, Clearly, P * 1 , which accounts for the potential energy in bonds and angles, is the AIC-best model. Note that this differs from the Bayesian implementation, in which P * 2 , consisting of bonds and Lennard-Jones 12-6 interactions, is chosen to be the most plausible.
In the validation scenario comprised of two chains of polyethylene, the MLE is updated, and we calculate using (23) and (24), respectively. Moving into the second validation scenario yields rendering the Category 2 model invalid, and necessitating another iteration through OPAL to Category 3. Note that when OPAL was implemented with Bayesian plausibilities as a model selection criterion, the Category 2 model was found to be not invalid.
Recalling that γ 1 measures only the difference in mean, while γ 2 takes into account the entire distribution of the potential energy, this result implies that the Bayesian approach better accounts for the various uncertainties present in this application. Continuing another iteration of OPAL, we select the Category 3 models such that M * = {P * 1 , P * 2 } = {P 11 , P 12 }. Both Category 3 models take into account bond, angle, and Lennard-Jones interactions; they differ only in the representation of the Lennard-Jones interaction. Calculating the Akaike weights, making the model with 9-6 Lennard-Jones interactions "better" than that with 12-6 Lennard-Jones interactions. In the first validation scenario, and in the second validation scenario, making the Category 3 model P * 2 not invalid for use in the prediction scenario. The potential energy distributions for all the "AIC-best" models in each of the three Occam Categories in the validation scenarios are given in Fig. 2.
A summary of the Bayesian implementation of OPAL results presented in [11] and the frequentist, AIC-based version of OPAL presented here is given in Table 2. Recall that the parameters were updated in the first validation scenario in both the deterministic and  Bayesian processes. The MLE calibration is much closer to the data than the parameter distributions produced by Bayesian calibration, as can be seen by a comparison of γ 1 and γ 2 for S v1 . Although most of the validation metrics computed in the first validation scenario for the MLE models are lower than those computed for the corresponding Bayesian models, there is a larger jump in these values as the complexity of the scenario increases (e.g., to the second validation scenario). Consider, for example, the validation metric values produced in Level 2. The relative change in γ 1 from S v1 to S v2 for the Bayesian plausibility is about 55%, while the change in AIC is about 20,000%. For γ 2 , this relative change is <1% for plausibility and nearly 8% for AIC. This may imply that the Bayesian models are more robust for extrapolation to more complex scenarios. It should be noted that these results depend on the data y that is used to calibrate the parameters, if Bayes' rule or maximum likelihood estimation is used. Theoretically, as the amount of data increases, the Bayesian posterior π(θ|d) and the MLE θ * of the truly best model converge to the true distribution or true value of the parameters, respectively [19][20][21]. It can be argued that, similarly, as the amount of data increases, Bayesian plausibilities and AIC values will converge to indicate the model that will best represent reality. Figure 3 provides plots of the dependence of the AIC values on the available data for models.

Concluding comments
On the basis of the sample calculations described in this work on the problem of validating coarse-grained models of atomistic systems, the Akaike and Bayesian criteria for model selection provide an efficient alternative to the more rigorous methods of Bayesian plausibility. Examples of implementations of the OPAL algorithm in model selection and validation suggest that the information-theoretic AIC selection procedures, as expected, seem to provide acceptable criteria for model selection. But in at least one case considered here, the best model selected by AIC differed from that pointed to by Bayesian plausibilities. From a practical point of view, even if the chosen model selection criteria fail to select the best model among a set of models proposed for a prediction, and if this model is invalid, this fact will be caught during the validation phase of OPAL. The better computational efficiency of the AIC methods in comparison with Bayesian plausibilities could make feasible new approaches to model selection and validation in the presence of uncertainties.
Author details 1 Sandia National Laboratories, Albuquerque, NM, USA, 2 Institute for Computational Engineering and Sciences, The University of Texas at Austin, Austin, TX, USA.