Evidential deep learning for interatomic potentials

Domain: physics
Published: 2025-12-20
Editorial outcome: major_revision

Abstract
Machine learning interatomic potentials have been widely used to facilitate large-scale molecular simulations with accuracy comparable to ab initio methods. To ensure the reliability of the simulation, the training dataset is iteratively expanded through active learning, where uncertainty serves as a critical indicator for identifying and collecting out-of-distribution data. However, existing uncertainty quantification methods tend to involve either expensive computations or compromise prediction accuracy. Here we show an evidential deep learning framework for interatomic potentials with a physics-inspired design. Our method provides uncertainty quantification without significant computational overhead or decreased prediction accuracy, consistently outperforming other methods across a variety of datasets. Furthermore, we demonstrate applications in exploring diverse atomic configurations, using examples including water and universal potentials. These results highlight the potential of our method as a robust and efficient alternative for uncertainty quantification in molecular simulations.

Introduction
Molecular dynamics (MD) simulation provides atomic insights into physical and chemical processes and has become an indispensable research tool in computational physical science1–3. Classical MD simulation uses an empirical potential function to determine interatomic forces4,5, which is computationally efficient but not accurate enough, especially when polarization or many-body interactions are important6. In contrast, the ab initio approach for modeling atomic interactions is based solely on fundamental physical principles, leading to generally higher accuracy and transferability7,8, but the high computational cost limits the size of systems that can be simulated. To achieve both efficiency and accuracy, machine learning interatomic potentials (MLIPs) have been proposed9–12, which allow to learn ab initio interatomic potentials and performing MD simulations with much lower computational cost. MLIPs have been successfully applied in the study of amorphous solid13, catalysis14, chemical reaction15, and more.

One of the primary challenges to MLIP-based MD simulations lies in the construction of the training dataset, which should include various configurations that may appear during the simulation. Inadequate training data will lead to decreased accuracy or even failure of the simulations16,17. This challenge limits the application of MLIP-based MD simulations. Active learning based on uncertainty quantification (UQ) plays a crucial role in constructing training sets for MLIPs18–21. During active learning, configurations with higher uncertainties are sampled to enrich the training set. This process usually needs to be repeated dozens or more times19, and the computational cost required for UQ could be considerable. Therefore, a robust yet efficient method for UQ is desired.

Results
Preliminary

Machine learning interatomic potential (MLIP). MLIPs are used to predict energy and forces within a given atomic configuration. For a system comprising N atoms, MLIPs typically take the atomic species Z 2ZN and coordinates R2RN×3 as input and output the total potential energy E. The forces F2RN×3 exerted on the atoms are derived by calculating the negative gradient of E with respect to the coordinates. The primary distinction among various MLIPs lies in the algorithm used to convert the input information into vectorized features that represent the local atomic environments. These features are designed to be invariant or equivariant under translation, rotation, and permutation.

Aleatoric and epistemic uncertainty. Two categories of uncertainty can be modeled in deep learning43. Aleatoric uncertainty arises from inherent noise in data labels. In the context of MLIPs, data labels are obtained from ab initio calculations. Although highly rigorous ab initio calculations could limit the aleatoric uncertainty, they are often computationally prohibitive in practice. As a result, MLIPs are often trained on multiple datasets of varying quality, which introduces label noise. For example, in the MPtrj dataset, aleatoric uncertainty may stem from inconsistent Hubbard U correction or varying convergence criteria26. In contrast, epistemic uncertainty is caused by a lack of knowledge, typically due to insufficient data. This type of uncertainty can be reduced by adding more training data, often through active learning. For the sake of simplicity, the term uncertainty in the following results refers to epistemic uncertainty, unless otherwise specified. We further discuss aleatoric uncertainty in Supplementary Note 2.

Evidential deep learning. Evidential deep learning is an efficient method to estimate the uncertainty of the results predicted by neural networks. Starting from a maximum likelihood perspective, the targets are assumed to be drawn from a Gaussian distribution but with unknown mean and variance (μ,σ2). A Gaussian prior is placed on the unknown mean μ and an Inverse-Gamma prior on the unknown variance σ2, leading to the Normal Inverse-Gamma distribution with a set of parameters m=(γ,ν,α,β). Neural networks are then trained to infer m, and the prediction, aleatoric, and epistemic uncertainty are calculated as35:

E[μ]=γ, E[σ2]= β/(α-1), Var[μ]= β/(ν(α-1)) (1)

Framework of eIP

As illustrated in Fig. 1, eIP consists of an MLIP block for predicting energy and force, and an evidential quantile regression block for UQ. The eIP framework extends a regular MLIP by feeding its equivariant features into a lightweight network that outputs prior parameters m. The loss function is central to eIP. We used a mean absolute error (MAE) loss for energy fitting, whereas force predictions are optimized by a composite loss that combines negative log-likelihood (NLL) with a regularization term. Further details are provided in the “Methods” section. In designing eIP, we considered the following points, which are indispensable to achieving robust performances.

Locality. In most MLIPs, the potential energy is calculated as the sum of atomic contributions, E= \sum_{i=1}^{N} E_i, with the model learning the mapping from the local environment of the atom i to E_i. Therefore, we estimate the uncertainty associated with E_i rather than the total potential energy E. However, we do not have the ground truth for E_i. Fortunately, we can adapt the atomic forces instead of E_i to estimate the uncertainty per atom.

Directionality. We attribute uncertainty in MLIP predictions to the inadequate learning of the local atomic environment. Consequently, this uncertainty is directionally dependent. Our model produces a separate uncertainty value for each Cartesian component of the atomic force, rather than a single value per atom or system. This is illustrated using a three-atom toy system in Supplementary Note 1, where the uncertainty varies with direction. Beyond MLIPs, the directional dependence is also crucial for predicting other non-scalar properties, such as dipole moments, dielectric tensors, and Hamiltonians. In the following experiments, we employ the equivariant backbone PaiNN44 to extract equivariant features and output the parameters of the Normal Inverse-Gamma prior. We also apply eIP to other equivariant backbones to demonstrate the generality, and the results are provided in Supplementary Note 5.

Quantile regression. The original evidential deep learning framework for regression35 assumes that the targets are drawn from a Gaussian distribution, which may not adequately describe the target distribution of MLIPs. To alleviate this limitation, we adopt the Bayesian quantile regression framework45 that combines evidential deep learning with quantile regression, enabling UQ without relying on Gaussian assumptions. Unlike traditional regression that minimizes mean squared error, our approach learns a specified quantile q of the force distribution. This is achieved by employing an asymmetric Laplace distribution, which is subsequently reformulated as a scale mixture of Gaussians to facilitate Bayesian inference. The computational procedure of Bayesian quantile regression is analogous to the original evidential framework. The key difference lies in the loss functions, which are specifically designed for quantile estimation. Our implementation details are provided in the “Methods” section.

Experiments

ISO17 dataset. We started by assessing the performance of eIP using the ISO17 dataset, which comprises MD trajectories of C7O2H10 isomers. This dataset is divided into in-distribution (ID) and out-of-distribution (OOD) subsets, making it particularly suitable for uncertainty quantification (UQ). In the ID scenario (known molecules/unknown conformations), the test molecules are also present in the training set. In contrast, the OOD scenario (unknown molecules/unknown conformations) involves test molecules that are not in the training set. The training set contains 400,000 conformations, which is a substantial amount for such small molecules. Therefore, we also explore the impact of training data volume. Specifically, we train the model using 1%, 5%, 30%, and 100% of the training data, respectively. Figure 2a–d shows the scatter plots that compare uncertainties with force errors for different amounts of training data, demonstrating positive correlations in both ID and OOD scenarios. The mean uncertainty and mean absolute error (MAE) for force predictions are shown in Fig. 2e, f, respectively. As expected, both metrics decrease with an increase in the amount of training data. Furthermore, we evaluated the reliability of UQ using additional metrics, including Spearman’s rank correlation coefficient and the area under the receiver operating characteristic curve (ROC-AUC), whose definitions are provided in the “Methods” section. As shown in Fig. 2g, h, both Spearman’s rank correlation coefficient and ROC-AUC improve as the amount of training data grows. In the ID scenario, Spearman’s rank correlation coefficients ranging from 0.74 to 0.86 and ROC-AUC values ranging from 0.86 to 0.93 indicate the strong performance of eIP. In the OOD scenario, although the test set molecules are absent from the training set, the evaluation metrics still provide reasonable uncertainty estimates.

Silica glass dataset. We then evaluate eIP’s performance for more complex systems using a silica glass dataset, which comprises large bulk structures. Given the challenges in partitioning large structures into ID and OOD datasets, we adopted the dataset partition scheme consistent with the previous study32. We also compare eIP with other UQ methods, including ensemble, Monte Carlo dropout, Gaussian mixture model (GMM), and Mean-variance estimation (MVE), whose implementations are provided in Supplementary Note 6. Figure 3a shows the scatter plots of uncertainties versus force errors and indicates that all methods achieve positive correlations. Figure 3b presents the computational efficiency analysis of the five methods. Despite the good performance of the ensemble method, it incurs higher computational costs during both training and inference due to the need for multiple independent MLIPs. Since our implementation uses four independent models, training and inference times are approximately four times longer than for a single model. It is worth noting that the actual time cost is highly dependent on implementation. In practice, the actual time required to train an ensemble can be less than the number of models implies. The Monte Carlo dropout method requires multiple inferences to obtain uncertainty, and the corresponding computational cost is similar to the ensemble method. GMM quantifies uncertainty through an iterative expectation-maximization algorithm and therefore requires additional computation time. Both MVE and eIP have minimal training and inference times, comparable to those of a normal MLIP. Regarding the force prediction accuracy shown in Fig. 3c, ensemble, GMM, and eIP achieve the lowest errors, while dropout and MVE exhibit larger errors. Figure 3d, e further illustrate the comparison of Spearman’s correlation and ROC-AUC, respectively. Notably, Fig. 3e shows that eIP performs even better than the ensemble method on the ROC-AUC metric.

Applications

Active learning with eIP. UQ plays a key role in active learning for training set construction. The quality of the training set is particularly crucial for MLIP, as the accuracy of MLIPs can significantly decrease when encountering unseen atomic configurations, leading to the collapse of simulations16. Figure 4a illustrates a typical active learning workflow for MLIPs, where the data points with high uncertainty are iteratively explored to enrich the training set. In addition, uncertainty-driven dynamics (UDD) simulation39 can be employed to enhance sampling efficiency. In UDD simulations, the potential energy surface is modified so that the atomic configurations with higher uncertainties are assigned lower potential energies, and consequently, these structures become more accessible, as indicated in Fig. 4b. The implementation of UDD simulation with eIP is provided in the “Methods” section.

We demonstrate the active learning process with eIP, using a water dataset as an example. In each iteration, we performed standard MD and UDD simulations using our eIP model. We refer to the UDD simulation driven by the uncertainty quantified by eIP as the eIP-UDD simulation. Changes in uncertainty and energy over simulation time are illustrated in Fig. 4c, where the uncertainty for each configuration is the mean atomic uncertainty determined by Eq. (7). The initial training set comprises 1000 configurations sampled from a classical MD simulation trajectory generated using an empirical force field. The abnormal energy fluctuations suggest that both the MD and eIP-UDD simulations collapse very early. After the first iteration, the MD simulation remains stable for the full 50 ps. Although the eIP-UDD simulation collapses after 20 ps, the uncertainty increases over time, indicating that more previously unseen configurations are explored during the eIP-UDD simulation. After the second iteration, both the MD and eIP-UDD simulations remain stable for the full 50 ps. We also observe that the uncertainty does not increase significantly, and this may suggest that configurations are explored sufficiently around certain local minima. The UDD parameters can be adjusted to further explore a broader configuration space. Details of the settings of UDD parameters are provided in Supplementary Note 7 and Supplementary Figs. 9–12.

Application of eIP in universal MLIP. Finally, we explored the performance of eIP in universal MLIPs. To this end, we trained the model on the Materials Project Trajectory (MPtrj) dataset46. The hexibin plots and the ROC curve in Fig. 5a–c demonstrate the performance of eIP on such a large dataset. We then conducted eIP-UDD simulations to test the performance of eIP in enhanced sampling. We selected two distinct materials as examples, namely lithium iron phosphate (LiFePO4) and polydimethylsiloxane (PDMS). LiFePO4 is a mature commercial cathode material for lithium-ion batteries, while PDMS is a widely applied organosilicon polymer material. These materials serve as benchmarks for evaluating the configurational sampling performance of eIP-UDD simulations for both inorganic crystalline and organic polymeric systems. To evaluate the diversity of the generated configurations, we calculated the configurational entropy for each trajectory, as detailed in the “Methods” section. For each material, changes in potential energy, uncertainty, and configurational entropy over simulation time are shown in Fig. 5d–i. In Fig. 5d, the initial LiFePO4 configuration was a structure from the Materials Project. No significant energy fluctuations were observed during the simulation. The brief rise in potential energy during the early stage of the eIP-UDD simulation occurs due to the modified potential energy surface that forces the system to escape the original local minimum. The PDMS configuration in Fig. 5g was built in an initial state, and the energy decrease corresponds to the structural relaxation process. In Fig. 5e, h, the trajectory of the eIP-UDD simulation has a higher uncertainty than that of the conventional MD simulation, as expected. The results of the configurational entropy in Fig. 5f, i further prove that the eIP-UDD simulations have obtained more diverse configurations.

UQ is a critical topic in various fields of machine learning, particularly in scientific applications such as molecular simulations based on MLIP. Conventional UQ methods suffer from either high computational costs or decreased prediction accuracy. In this work, we propose a single-model UQ method, called eIP, which achieves both efficiency and accuracy, as demonstrated by extensive experiments in various applications. The eIP framework incorporates locality, directionality, and quantile regression, all of which are essential for achieving optimal results. This is evident from the ablation study presented in Supplementary Note 3, where the absence of any single component leads to a noticeable decline in performance.

Although ensemble methods have been widely used in active learning, they typically require training four or more models simultaneously. In practice, this process usually involves dozens or more iterations and takes a significant amount of time and computational resources to obtain a satisfactory training set. As a result, single-model UQ methods, such as eIP, have the potential to save several months in applications, making eIP a more efficient alternative when time constraints and computational resources are a significant concern. In addition, for large-scale simulations, ensemble methods require a significant amount of computation to evaluate the reliability of MLIP-based MD simulations, while eIP facilitates real-time assessment without incurring noticeable additional costs.

Methods
Formalism of eIP

We employ quantile regression with maximum likelihood estimation to better model the uncertainty of MLIPs. Quantile regression is solved by minimizing the tiled loss for a given quantile q:

L_q(ϵ) = max(qϵ, (q-1)ϵ), (2)

where ϵ is the error. The quantile loss is equivalent to a maximum likelihood estimation for an asymmetric Laplace distribution with mean μ, variance σ, and an asymmetrical parameter equal to the quantile q47. The likelihood function can be expressed as a scalar mixture of Gaussians48,49 N(μ+τz, ωσz), where τ = (1-2q)/(q(1-q)), ω = 2/(q(1-q)), z = exp(1). We assume that the atomic forces F2RN×3 come from a Gaussian distribution, but the mean and variance are unknown. For instance, the x-component of the force on the atom i follows:

f_ix ~ N(μ_ix+τz_ix, ωσ_ixz_ix) (3)

By placing a Normal-Inverse-Gamma (NIG) evidential prior p(μ_ix, σ^2_ix | m_ix) with a set of parameters m_ix=(γ_ix, ν_ix, α_ix, β_ix)35,45. As a result, γ_ix is equal to the predicted force

E[μ_ix] = \int \int μ_ix p(μ_ix) dμ_ix = γ_ix, (4)

and the x-component of epistemic uncertainty for the atom i is

Var[μ_ix] = \int \int μ^2_ix p(μ_ix) dμ_ix - (E[μ_ix])^2 = β_ix / (ν_ix(α_ix-1)) (5)

The y- and z-components are computed similarly. We define the uncertainty σ_i associated with the atom i as

σ_i = \sqrt{β_ix / (ν_ix(α_ix-1)) + β_iy / (ν_iy(α_iy-1)) + β_iz / (ν_iz(α_iz-1))} (6)

The uncertainty for a configuration composed of N atoms is determined by computing the average:

σ = 1/N \sum_{i=1}^{N} σ_i (7)

It should be noted that averaging can lead to a loss of local information. While the maximum uncertainty value is an alternative, it is susceptible to intrinsic errors in uncertainty quantification. The weighted quantile is a more robust choice, but it requires careful selection of the weights.

The model is trained by maximizing the probability p(f_ix | m_ix). Marginalizing out the likelihood parameters μ_ix and σ^2_ix:

p(f_ix | m_ix) = \int_{σ^2_ix=0}^{\infty} \int_{μ_ix=-\infty}^{\infty} p(f_ix | μ_ix, σ^2_ix) p(μ_ix, σ^2_ix | m_ix) dμ_ix dσ^2_ix (8)

By placing the NIG prior on μ_ix and σ^2_ix, this integral has an analytical solution:

p(f_ix | m_ix) = St(f_ix; γ_ix+τz_ix, 2β_ix(1+ν_ix+α_ixωz_ix)/(ν_ixα_ix), 2α_ix) (9)

where St(f; μ_St, σ^2_St, v_St) is the Student t-distribution evaluated at f with location parameter μ_St, scale parameter σ^2_St, and v_St degree of freedom. We then obtain the negative log-likelihood (NLL) loss function45:

L_NLL = -log p(f_true_ix | m_ix) = 1/2 log(π/(ν_ixα_ix)) - 1/2 log(Γ(α_ix)/Γ(α_ix+1/2)) + (α_ix+1/2) log(Ω(f_true_ix - (γ_ix+τz_ix))^2 + 2ν_ix + Ω) + log Γ(α_ix) (10)

where Ω = 4β_ix(1+ωz_ixν_ix), z_ix = α_ix/β_ix, and Γ(⋅) is the gamma function. We use an evidence regularizer35 so that the model tends to output low confidence when the predictions are incorrect:

L_R = ρ_q(f_true_ix - γ_ix) ⋅ Φ_ix (11)

where Φ_ix = 2ν_ix + α_ix + 1. The model learns to reduce confidence by outputting lower values for ν and α, or a higher value for β. Consequently, as demonstrated by the quantitative results in Supplementary Table 1, this regularization term effectively mitigates overconfidence. The y- and z-components are computed similarly. Finally, the overall loss function, including the L1 loss for energy prediction, is:

L = |E_true - E| + w \sum_{i=1}^{N} \sum_{a \in {x,y,z}} L_NLL + λL_R (12)

where w and λ are hyperparameter to adjust the weighting of each term. The effects of these parameters on the results are discussed in Supplementary Note 4.

The eIP model in this work is implemented on the PaiNN44 backbone as an example, but it is also applicable to other equivariant backbones (see Supplementary Note 5). In contrast to the standard PaiNN model, the eIP model incorporates an additional evidential block, which takes the equivariant features as input and produces the output α, β, and ν, as illustrated in Supplementary Fig. 7. Since the evidential block is lightweight compared to the message-passing layers in the backbone, the additional computational overhead is minimal. The model parameters are trained by minimizing the overall loss function Eq. (12).

Datasets

ISO17 dataset. The ISO17 dataset50 was obtained from http://quantum-machine.org/datasets/. We adopted the original splitting strategy for the training, validation, and test sets. For training sets of different sizes, the smaller training sets were randomly sampled from the largest training set, containing 400,000 conformations.

Silica glass dataset. The silica glass dataset is obtained from a previously published study32. It comprises 1691 configurations, each containing 699 atoms (233 Si and 466 O atoms). These configurations were generated through MD simulations with a force-matching potential51 under various conditions, followed by density functional theory (DFT) calculations to obtain energies and forces. We adopted the same dataset splitting scheme as described in the ref.32. Partitioning these structures into ID and OOD datasets is challenging, as it is difficult to find configurations with atomic environments entirely distinct from one another. To reflect a more generalized evaluation under more extreme conditions, the training set includes only structures generated under low-temperature, low-deformation-rate conditions, while the test set contains structures extracted from trajectories at higher temperatures and higher deformation rates.

Water dataset. The initial water training set is taken from our previous work17. It comprises 1000 configurations sampled from classical MD trajectories with the SPC/E force field52. Each configuration contains 288 atoms with periodic boundary conditions. During active learning, we ran UDD simulations at 300 K and sampled 1000 configurations for each iteration. The energies and forces are determined using density functional theory (DFT) calculations employing the cp2k software package53 with the PBE-PAW-DFT-D3 method54–56.

MPtrj dataset. The MPtrj dataset46 is a collection of MD trajectories designed for training a universal potential. It comprises millions of configurations covering 89 elements, and the energies and forces are determined using DFT calculations. We adopted the original splitting strategy with an 8:1:1 training, validation, and test ratio.

Evaluation metrics

Spearman’s rank correlation coefficient. Spearman’s rank correlation is a non-parametric measure of the strength and direction of association between two ranked variables. Unlike Pearson’s correlation, which accesses linear relationships, Spearman’s rank correlation evaluates how well the relationship between two variables can be described using a monotonic function. We expect a larger error to be associated with a higher uncertainty, and their correlation does not necessarily need to be linear. Therefore, Spearman’s rank correlation coefficient was used to assess the reliability of the uncertainty. A coefficient of 1 means perfect correlation, and a coefficient of 0 indicates that there is no correlation between the ranks of the two variables.

Area under the receiver operating characteristic curve. The receiver operating characteristic (ROC) curve is a graphical representation of a classifier’s performance. The area under the ROC curve (ROC-AUC) provides a complementary evaluation metric for UQ that avoids the possible limitations of using Spearman’s rank correlation coefficient alone. Following the approach of a previous study32, we designed a classification task in which predictions with high errors are expected to exhibit high levels of uncertainty. An error threshold (ε_c) and an uncertainty threshold (U_c) are defined to classify data points. A data point is classified as a true positive (TP) if both its true error and estimated uncertainty exceed their respective thresholds (ε > ε_c and U > U_c); a false positive (FP) if the error is below its threshold but the uncertainty is above (ε ≤ ε_c and U > U_c); a true negative (TN) if both are below their thresholds (ε ≤ ε_c and U ≤ U_c); and a false negative (FN) if the error is above its threshold but the uncertainty is below (ε > ε_c and U ≤ U_c). We set the threshold to be at the 20th percentile as in ref.32. The ROC-AUC score ranges from 0 to 1, with a score of 1 denoting a perfect classifier and 0.5 indicating performance no better than random choice.

Configurational entropy. Configurational entropy57,58 quantifies the number of ways that atoms in a system can be arranged. High entropy indicates that the system is likely to take on many different arrangements, whereas low entropy implies a more ordered, less random state. We used configurational entropy as a metric to measure the diversity of configurations obtained during MD and UDD simulations. The configurational entropy is defined as the Shannon entropy59,60:

S_conf = - \sum_{i} p_i log p_i, (13)

where p_i is the probability of the system being in state i. We then projected states onto a discretized order parameter grid and calculated the frequency of these order parameters within a simulation trajectory. For LiFePO4, the selected order parameters were the P-O-Fe angle and the PO4 tetrahedral distortion. For PDMS, we selected the end-to-end distance and the radius of gyration as the order parameters. To determine the probability distribution, the order parameter space was discretized into an N_e × N_e grid, and the frequency of configurations within each grid cell was calculated. The configurational entropy was normalized by dividing it by the maximum possible entropy value, 2 log(N_e), resulting in values between 0 and 1. A larger grid size N_e offers a finer resolution but may suffer from statistical noise, while a smaller N_e provides more robust statistics at a lower resolution. We used N_e = 40 for all reported results. Varying the value of N_e does not significantly affect the results, as the configurational space was sampled sufficiently in our simulations.

Molecular dynamics (MD) simulations

MD simulations were performed using the Atomic Simulations Environment (ASE) Python library61. The simulations are set with a timestep of 0.1 fs in the canonical (NVT) ensemble. The Berendsen thermostat62 was used with a coupling temperature of 300 K and a decaying time constant τ of 100 fs. The atomic velocities were initialized according to the Boltzmann distribution at 300 K. The initial water configuration was selected from the water test set. The LiFePO4 configuration was obtained from the Materials Project, comprising 168 atoms in the unit cell. The PDMS configuration was constructed using three polymer chains with a polymerization degree of 25 and a density of 0.97 g⋅cm−3, containing 759 atoms in total. All systems were modeled with periodic boundary conditions.

Uncertainty-driven dynamics (UDD) simulations

The UDD simulation technique utilizes a bias energy that favors configurations with higher uncertainties. Kulichenko et al. introduce a bias energy39 defined as:

E_bias(σ^2) = A exp(-σ^2/B) - 1, (14)

where the parameters A and B are chosen empirically. The bias force F_bias is then determined by calculating the negative gradient of the bias energy:

F_bias = -∇(E_bias(σ^2)) = -E_bias(σ^2) * ∇σ^2, (15)

By leveraging eIP for UQ, the gradient of σ can be obtained through automatic differentiation.

Notably, the bias force could become exceptionally large, leading to the collapse of molecular simulations. We found that limiting the magnitude of the bias forces using a clipping strategy proved not effective. To prevent this issue, we incorporate a Gaussian term to limit the magnitude of the bias force with two additional empirically chosen parameters C and D:

F_limited_bias = F_bias * exp(-F_bias^2 / (2πC * 2C^2)) (16)

This adjustment of bias force implies a new bias energy formulation and ensures more stable UDD simulations. Detailed discussions about the empirical parameters A, B, C, and D are provided in the Supplementary Note 6. Finally, the combined force F + F_limited_bias is used to guide the simulations toward configurations with higher uncertainties, enhancing the sampling for more diverse atomic configurations.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The manuscript presents a novel approach to evidential regression for machine-learned interatomic potentials. A key\nstrength of the proposed method is that it enables single-model uncertainty quantification (UQ), making it computationally\nefficient and well-suited for large-scale applications. The authors adapt the output of the PaiNN architecture to their\nevidential UQ framework, allowing it to predict an uncertainty value for each force component. Notably, the model’s\nepistemic uncertainty correlates well with the actual force prediction error.\nThe authors demonstrate the effectiveness of their evidential interatomic potential (eIP) on standard molecular benchmark\nsystems, showing that it performs on par with, or better than, other single-model UQ approaches. In particular, the method\nshows improved performance over Monte Carlo dropout and achieves comparable accuracy to computationally more\ndemanding ensemble-based UQ.\nThe manuscript is well written and represents a meaningful contribution to the field of uncertainty quantification for machine-\nlearned force fields. To strengthen the work and broaden its impact, I encourage the authors to consider the following main\npoints:\n1. Comparison to recent UQ methods: The benchmarking shown in Figure 3 appears to be limited to the methods evaluated\nin earlier work (primarily Ref. 30). However, several newer approaches to single-model UQ have recently been introduced,\nincluding LTAU (https://doi.org/10.1088/2632-2153/adb4b9), Prediction Rigidity (https://doi.org/10.1088/2632-2153/ad805f),\nand Shallow Ensembles (https://doi.org/10.1088/2632-2153/ad594a). Including a comparison to these more recent methods\n—or at least discussing them—would help position the present work with respect to current work and highlight its relative\nstrengths.\n2. Architecture transparency and model generality: A schematic of the uncertainty output layer (perhaps included in the\nSupporting Information) would help make the model more accessible and easier to implement for other groups. Additionally,\nit would be informative to test the generality of the eIP framework by applying it to other modern equivariant neural network\narchitectures beyond PaiNN. Even a small-scale experiment comparing eIP to ensembling across different architectures\nwould help support the claim that this approach is robust and generally applicable.\n3. Focus of the active learning experiments: The section comparing different active learning strategies feels somewhat\ndisconnected from the core contribution of the paper. Rather than testing multiple acquisition schemes, it might be more\ninsightful to fix a single active learning strategy and compare how different UQ methods—e.g., eIP, ensembles, MC dropout\n—guide data acquisition and model improvement. This would more directly demonstrate the utility of the new UQ approach,\nsimilar to the strategy in Ref. 30. As it stands, the fact that eIP-UDD does not clearly outperform vanilla MD sampling may\nreflect limitations of the UDD scheme rather than the UQ model itself.\nAdditional comments, ordered by appearance in the manuscript:\nLines 91–92: The term uncertainty-driven dynamics (UDD) is introduced but not sufficiently explained. Moreover, the\nphrasing “active learning and UDD” is somewhat confusing, as UDD is generally considered a form of active learning. It\nwould be helpful to clarify the distinction the authors intend here. The authors might also consider citing additional recent\nUDD variants, such as: https://doi.org/10.1038/s41524-023-01104-6, https://doi.org/10.1038/s41524-024-01254-1,\nhttps://doi.org/10.1063/5.0246178\nLines 93–94: The statement regarding a “universal potential that achieves real-time UQ” could be clarified. Do the authors\nmean that UQ is possible concurrently with simulation? If so, this is not unique to the present method—techniques such as\nensembling, MVE, and GMM can also provide real-time uncertainty estimates.\nLine 100: Should the atomic species be denoted as \"Z in Z^N\" rather than \"Z in R^N\"?\nLines 136–137: The explanation of “directionality dependence” is too brief to be easily understood without referring to\nSection S1. If I understand correctly, the method produces not a single uncertainty per atom or system, but one per Cartesian\ncomponent of the atomic force. This is a key feature and should be explained more clearly in the main text. The authors\nshould also highlight the significance of this directionality—for example, its novelty or potential applications.\nFigure 1: While the structure of the loss function is described in the Methods section, a short explanation in the main text\nwould help readers appreciate the role of each term and the conceptual difference from standard regression approaches.\nAdditionally, the notation in Figure 1 (using theta) appears inconsistent with Equation 10 (which uses \"w\"), assuming that\nw=1/θ. This should be clarified.\nLines 144–146: While most readers will be familiar with Gaussian distributions and conventional UQ, the mention of\nBayesian quantile regression may need elaboration. A brief explanation of what this model entails and its relevance to the\nauthors' approach would improve accessibility.\nLine 172 and Lines 323–325: The term “favorable range” is used without definition. Could the authors specify what is meant\nby a favorable versus unfavorable AUC range? Additionally, to help readers interpret Figure 2h, it would be helpful to define\nin the main text—or clearly reference the relevant part of the Methods section—how a classification test is constructed from a\nregression task. This is key for understanding the figure without needing to study Ref. 30. Specifically, could the authors\nbriefly explain how true/false positive/negative cases are determined?\nLine 182: The phrase “four times the training time” is misleading without context. This scaling is due to the use of an\nensemble of size 4, which does not appear to be mentioned earlier. It would be more accurate to state that training time\nscales linearly with the number of models in the ensemble. Similarly, inference time for both ensembles and MC dropout (a\npseudo-ensemble) also scales linearly with the ensemble size—here, 4.\nLine 186: The authors mention that expectation-maximization (EM) increases inference cost for the GMM model. However,\nEM is typically used during GMM training, not inference. At prediction time, GMM inference should be relatively lightweight.\nCould the authors clarify what is contributing to the increased cost of GMM inference, and why it is reported as slower than\nan ensemble?\nFigure 3c: Could the authors consider reporting the force RMSE instead of the MAE? For vector quantities, the L1 norm is\nnot invariant under rotation, so MAE can yield inconsistent values depending on molecular orientation. The RMSE (L2 norm)\ndoes not suffer from this problem and would be more appropriate here.\nFigure 4c: Up to this point, the authors emphasize that their model yields directional uncertainties. However, in Figure 4c, the\ny-axis is labeled simply as uncertainty. Since the model appears to predict a tensorial uncertainty per configuration, it would\nbe helpful to explain how this scalar value is derived from the full prediction. A brief explanation or a pointer to the relevant\nequation in the Methods section (where it is explained) would clarify this.\nLines 200–203: This brief explanation of how directional uncertainty helps identify chemical reactivity is insightful. Consider\nmoving a version of this explanation to the Introduction, where it could better motivate the relevance of the method.\nLines 206–207: The phrasing “standard MD and eIP-UDD” is potentially confusing. Since both simulation types are carried\nout using the trained eIP model, I suggest rephrasing to something like “we performed both standard MD and UDD\nsimulations using our eIP model”. Alternatively, if the term eIP-UDD is retained throughout the manuscript, it should be\nclearly defined at first occurrence to avoid the impression that two different force fields are being compared.\nLine 215: The statement that the simulations are “stable after 50 ps” might be misleading. Do the authors mean that the\nsimulations remain stable for the full 50 ps? Also, the current phrasing does not communicate what advantage the UDD\nsimulation provides, especially since both UDD and standard MD lead to stable models after two active learning iterations.\nCould the authors clarify the added value of UDD in this context?\nFigure 5: The caption should specify which data set is used. In panels (d) and (g), the abrupt changes in potential energy at\nthe start of the trajectories are concerning. However, they could plausibly be attributed to the lack of proper thermalization.\nWas there a heating or equilibration stage prior to the MD runs shown here? If not, a brief explanation of why this was\nomitted—or a clarification in the figure caption—would help interpret the results.\nLines 228 and 333–334 (configurational entropy): A reference to the relevant Methods or SI section would be helpful here.\nAlso, is there a standard citation for the specific approach used? Furthermore, could the authors expand on how\nconfigurational entropy is computed in this work? For example, is a long unbiased MD simulation first performed to collect\nhistograms of key order parameters, which are then used to estimate configuration probabilities? Or are the probabilities\ninferred from the training set? If so, this method relies on the assumption of a sufficiently unbiased and ergodic sampling,\nwhich may not hold in the case of data sets obtained through active learning. In such scenarios, the resulting histograms\nmay not reflect a proper Boltzmann distribution, making entropy estimates difficult to interpret.\nLine 258: When referring to \"residue\", do the authors mean \"residual\" as in the difference between the predicted and\nreference values?\nEquation 7: Have the authors considered using alternative aggregation functions for estimating the global uncertainty,\nbeyond the mean? For instance, in larger systems like silica, it is possible that a few atoms experience large force prediction\nerrors, which could lead to failure in the simulation—even though the average uncertainty remains low. A maximum or\nweighted quantile might better reflect risk in such edge cases.\nAlso, should both sides of the equation include a square term? As written, it appears the square of the per-atom uncertainties\nmay be missing.\nEquations 8 and 9: The motivation and derivation for these equations could be expanded. As it stands, they are difficult to\nfollow for readers unfamiliar with Bayesian quantile regression. A brief explanation of how these formulations emerge from\nthe underlying assumptions or loss functions would improve accessibility.\nLines 250 / 351: Is it advisable to start simulations directly at 300 K, particularly in early active learning iterations when the\nMLIP may still be unreliable? This practice could increase the likelihood of simulation instabilities or system blow-ups.\nCould the authors comment on this choice and whether a lower initial temperature or a temperature ramp was considered?\nVery minor recommendation regarding SI front page: It would be helpful to include a short title page for the SI, stating the title\nof the paper, author names, and affiliations. This is standard in many journals and aids in readability and citation.\nSection S3: The sentence “When ablating locality, we chose global embedding to estimate the uncertainty associated with\nthe energy.” is unclear. There is no “global embedding” in the PaiNN architecture. If the authors refer to pooling features\nacross atoms or to removing the spatial locality constraint in uncertainty prediction, they should reformulate this more\nprecisely.\nSection S4: “In contrast to the standard PaiNN model, the eIP model incorporates additional fully connected layers.” It is\nunclear whether the authors refer to changes in the output head or in the message passing layers themselves. The original\nPaiNN already defines multiple output architectures, such as simple feed-forward networks for scalar outputs and gated\nequivariant blocks for vector-valued quantities. If the authors modified or replaced these, it would be valuable to include a\nvisual diagram of the architecture, especially to clarify whether these changes are post-message-passing or occur in\nparallel.\nSI Eq. 5 and Eq. 7: Could the authors specify the dimensionality or size of the term ε?\nSI Eq. 9 vs Main Text Eq. 6: There appears to be an inconsistency between Eq. 9 in the SI and Eq. 6 in the main text. In the\nSI, the variances are summed directly, while in the main text, the squares of the variances are added before taking the\nsquare root. The authors should clarify this discrepancy and ensure consistency between formulations.\nSI, page 9 (\"final layer\"): The term “final layer” could refer to the last layer of the message passing network or to the last layer\nof the output head. Please clarify which is meant here, especially as this distinction may affect how uncertainty is computed\nor aggregated.", "reviewer_2": "This report details the use of evidential deep learning to capture epistemic uncertainty in MLIPs at reduced computational\ncost as compared to standard methods such as ensembling, Monte Carlo dropout, GMM, and MVE. Uncertainty-driven dynamics\nsimulations are used alongside their eIP method in an active learning process. Their eIP method is of significance\nto the field of MLIPs and will help to model epistemic uncertainties without the need for expensive model ensembling.\nIn the discussion of aleatoric versus epistemic uncertainty, while it is true that strict curation of a dataset would limit the\naleatoric uncertainty, this type of careful data curation is often not what is done in practice, especially in newer foundation\nMLIPs that often combine multiple large datasets. It should be noted in the main text that aleatoric uncertainty is present in\nmany pretrained foundation/universal models. In particular, this should be discussed for the MPtrj dataset, as there are\nknown sources of aleatoric uncertainty in that dataset. See this report for example: https://www.nature.com/articles/s41524-\n025-01572-y\nModel ensembling is known to produce overconfident uncertainties during inference. The authors use an evidence\nregularization function to reduce model overconfidence. More discussion should be given concerning this topic. A\ndemonstration should be included to show the effect of the evidence regularization function on reducing overconfident\nuncertainties.\nBecause the silica glass dataset used in this work is not commonly used in benchmarking MLIPs, more information is\nneeded describing the dataset, including the differences between the ID and OOD datasets partitioned in reference 30.\nPlease provide the empirical force field used to generate the classical MD water dataset used during eIP-UDD.", "reviewer_3": "Review for “Evidential Deep Learning for Interatomic Potentials”\nIn the manuscript “Evidential Deep Learning for Interatomic Potentials” the authors Han Xu et al describe a method for\ntraining machine learning interatomic potentials to predict parameters in a distribution rather than a singular energy and\nforce value. This particular parameters inferenced in this distribution are then used to indicate whether a specific atomic\nconfiguration in active learning should be simulated with electronic structure theory or not. The authors due a commendable\njob of comparing this method to other uncertainty quantification methods and link it to previously developed accelerated\nsampling techniques (UDD). In general, I find the manuscript to be well can clearly written and provide a valuable\nmethodology for reducing the computational cost associated with active learning. My only concern is that the technique\ndemonstrated here is only of use to those in the community actively building MLIP training data with UQ methods. This is a\nrelatively small community and thus this paper may not be suitable for nature communications simply because of the niche it\noccupies. This is a decision to be made by the editor.\nMy only real objection in the manuscript is that the ensemble method is shown to be a clear factor of 4 more expensive for\ntraining and inference tasks. This is not quite correct for various practical and technical reasons. For training, most HPC\nnodes have multiple GPUs. Due to the limitation of data transmission/memory synchronization, training 4 separate models\non 4 separate GPUs only takes about 2 times longer than training 1 model on 4 GPUs. This is obviously package/model\ndependent but would guess similar trends would hold for other Pytorch based models. Further, in inference, ensemble-\nbased methods, when properly implemented, are not “# of networks” more expensive than a single model. This is because\nthings like mapping of atoms to symmetry functions (which are shared between models) can be re-used for each model. The\nrelatively cheap tensor operations of the MLIP are the only things that need to be repeated for each element of the ensemble.\nThis practically saves 40% of calculation time in the ensemble, making it much more competitive than Fig 3 suggests. I\nrealize these details are difficult to discuss in a general sense, but possibly a footnote on ensemble timings is worth\nincluding.\nThis manuscript also suggests interesting questions outside of the present scope: Does aleatoric uncertainty actually relate\nto noise in the data labels? This is easy to explore in an electronic structure context because the convergence threshold of\nbasically every calculation is adjustable. It would be interesting to numerically verify that the author’s label of aleatoric\nuncertainty goes up when lower convergence thresholds are used. Further, does the ability to train a specific aleatoric\nuncertainty parameter enable the model to more efficiently use poorly converged simulations, potentially saving time on\nelectronic structure calculations?"}, "author_response": "Thank you for highlighting these recent works. We have carefully studied the mentioned\nmethods, namely LTAU-FF, LLPR (prediction rigidity), and DPOSE (shallow ensembles).\nWhile they offer novel and valuable solutions for uncertainty quantification (UQ), our eIP\nframework provides distinct advantages in theoretical foundation and practicality:\n- LTAU-FF is a very inspiring method that efficiently extracts uncertainty from training error\ntrajectories. However, its strong dependence on training dynamics can risk biased\nuncertainty estimates. We also noted that its reported performance decreases compared\nto the ensemble method, particularly on out-of-distribution (OOD) data.\n- LLPR quantifies uncertainty by perturbing model predictions and measuring loss\nsensitivity within a complete mathematical framework. Nonetheless, its practical\napplication relies on several non-trivial approximations, such as the use of generalized\nGauss-Newton (GGN) Hessian to avoid calculating second-order derivatives, the finite-\nwidth approximation in neural network Gaussian processes (NNGPs) and neural tangent\nkernel (NTK) theory, and the last-layer approximation. While this method has proven\neffective for predicting potential energy on benchmarks like QM9, its validity for MLIPs\n(especially for large models with millions of parameters) requires further verification on\natomic forces, which is more challenging due to their nature as energy derivatives.\n- DPOSE combines the ideas of ensemble and mean-variance estimation (MVE). Weight-\nsharing is used to reduce computational cost, and this is similar to the last-layer\napproximation in LLPR. However, this approximation may worsen the known issue of\nensemble overconfidence, especially in larger models. This concern is supported by\nexisting research (https://www.nature.com/articles/s41524-025-01572-y), which observed\nnotable epistemic overconfidence in a foundation model using a similar approach.\nDPOSE also incorporates techniques such as post-calibration and the continuous ranked\nprobability score (CRPS) to improve performance. Furthermore, it introduces an\nuncertainty propagation algorithm that estimates the uncertainty in derived properties.\nThis feature is particularly valuable for analyzing the results of molecular dynamics\nsimulations performed with MLIPs. We have also implemented DPOSE and compared it\nwith ensemble and eIP.\nIn contrast, our eIP framework provides an end-to-end solution for quantifying uncertainty in\nvectorized targets such as atomic forces. It analytically derives both aleatoric and epistemic\nuncertainties from an evidential prior, whose parameters are predicted by a lightweight\nnetwork. Crucially, eIP's performance is robust to model size and training data volume,\nmaking its advantages particularly significant for modern foundation models. The revised\nintroduction (lines 70-80) includes discussions of these recent methods. A comparison with\nthe shallow ensemble method is provided in the revised Supplementary Section S5.\n2. Architecture transparency and model generality: A schematic of the uncertainty output\nlayer (perhaps included in the Supporting Information) would help make the model\nmore accessible and easier to implement for other groups. Additionally, it would be\ninformative to test the generality of the eIP framework by applying it to other modern\nequivariant neural network architectures beyond PaiNN. Even a small-scale\nexperiment comparing eIP to ensembling across different architectures would help\nsupport the claim that this approach is robust and generally applicable.\nResponse:\nThank you for the suggestion. We have added a figure illustrating the architecture of the eIP\noutput layer (Figure S7). We also agree that demonstrating the generality of the eIP\nframework beyond PaiNN is crucial. In response, we applied the eIP framework to two recent\narchitectures: TorchMD-Net and AlphaNet. We evaluated their performance on the silica\ndataset and made comparisons between eIP, ensemble, and shallow ensemble (DPOSE).\nThe results shown in the revised Supplementary Section S5 strongly support our claim\nregarding the general applicability of eIP. The code has been uploaded to figshare\n(https://figshare.com/articles/dataset/Evidential_Deep_Learning_for_Interatomic_Potential/28\n805819)\n3. Focus of the active learning experiments: The section comparing different active\nlearning strategies feels somewhat disconnected from the core contribution of the\npaper. Rather than testing multiple acquisition schemes, it might be more insightful to\nfix a single active learning strategy and compare how different UQ methods—e.g.,\neIP, ensembles, MC dropout—guide data acquisition and model improvement. This\nwould more directly demonstrate the utility of the new UQ approach, similar to the\nstrategy in Ref. 30. As it stands, the fact that eIP-UDD does not clearly outperform\nvanilla MD sampling may reflect limitations of the UDD scheme rather than the UQ\nmodel itself.\nResponse:\nThe primary goal of this experiment was to demonstrate the application of eIP in active\nlearning, with a focus on showing that eIP-UDD can indeed sample more high-uncertainty\nconfigurations. As observed from the simulation trajectories, the eIP-UDD simulations\nsuccessfully explore configurations with higher uncertainty. Although the eIP-UDD simulations\nbecome more prone to instability, the results after the second iteration show that the\nsimulation remains stable. At this stage, the configurational space can be further expanded by\nadjusting simulation conditions (e.g., temperature, pressure) or fine-tuning UDD parameters\n(see revised Supplementary Section S7 and Figures S9-S12). We have revised the analysis\nof Figure 4 in the main text (lines 245-254).\nRegarding the comparison of eIP with other methods, a comprehensive evaluation has been\nconducted using the silica dataset (Figure 3). We note that directly comparing the\nperformance of various UQ methods within an active learning setting is challenging. Factors\nsuch as the splitting of training/test sets and various randomness in MD simulations can\nsignificantly influence the results. When the initial training set is small, newly added data\ngenerally improve model performance regardless of the UQ method used, and thus the\ndifferences between UQ methods are often minimal, just as shown in Ref. 30. We believe that\nthe choice of UQ method becomes more critical when the training set is already large. An\ninaccurate UQ method may fail to identify truly high-uncertainty data, leading to the selection\nof redundant data points.\nAdditional comments, ordered by appearance in the manuscript:\n4. Lines 91–92: The term uncertainty-driven dynamics (UDD) is introduced but not\nsufficiently explained. Moreover, the phrasing “active learning and UDD” is somewhat\nconfusing, as UDD is generally considered a form of active learning. It would be\nhelpful to clarify the distinction the authors intend here. The authors might also\nconsider citing additional recent UD variants, such as:\nhttps://doi.org/10.1038/s41524-023-01104-6, https://doi.org/10.1038/s41524-024-\n01254-1, https://doi.org/10.1063/5.0246178\nResponse:\nThank you. We have revised the manuscript accordingly (lines 99-105).\n5. Lines 93–94: The statement regarding a “universal potential that achieves real-time\nUQ” could be clarified. Do the authors mean that UQ is possible concurrently with\nsimulation? If so, this is not unique to the present method—techniques such as\nensembling, MVE, and GMM can also provide real-time uncertainty estimates.\nResponse:\nThank you. We demonstrate eIP's ability for concurrent UQ during simulations. While other\ntechniques can also achieve this, eIP offers advantages in efficiency and reliability. We have\nrevised the manuscript accordingly (lines 99-105).\n6. Line 100: Should the atomic species be denoted as \"Z in Z^N\" rather than \"Z in\nR^N\"?\nResponse:\nThank you. This has been corrected.\n7. Lines 136–137: The explanation of “directionality dependence” is too brief to be easily\nunderstood without referring to Section S1. If I understand correctly, the method\nproduces not a single uncertainty per atom or system, but one per Cartesian\ncomponent of the atomic force. This is a key feature and should be explained more\nclearly in the main text. The authors should also highlight the significance of this\ndirectionality—for example, its novelty or potential applications.\nResponse:\nThank you. We have revised the manuscript accordingly (lines 153-163).\n8. Figure 1: While the structure of the loss function is described in the Methods section,\na short explanation in the main text would help readers appreciate the role of each\nterm and the conceptual difference from standard regression approaches.\nAdditionally, the notation in Figure 1 (using theta) appears inconsistent with Equation\n10 (which uses \"w\"), assuming that w=1/θ. This should be clarified.\nResponse:\nThank you. We have revised the manuscript accordingly (lines 139-146). The notation \"theta\"\nin Figure 1 is replaced with \"w\" for consistency.\n9. Lines 144–146: While most readers will be familiar with Gaussian distributions and\nconventional UQ, the mention of Bayesian quantile regression may need elaboration.\nA brief explanation of what this model entails and its relevance to the authors'\napproach would improve accessibility.\nResponse:\nThank you. We have revised the manuscript accordingly (lines 164-175).\n10. Line 172 and Lines 323–325: The term “favorable range” is used without definition.\nCould the authors specify what is meant by a favorable versus unfavorable AUC\nrange? Additionally, to help readers interpret Figure 2h, it would be helpful to define in\nthe main text—or clearly reference the relevant part of the Methods section—how a\nclassification test is constructed from a regression task. This is key for understanding\nthe figure without needing to study Ref. 30. Specifically, could the authors briefly\nexplain how true/false positive/negative cases are determined?\nResponse:\nWe acknowledge that the word \"favorable\" is not quite appropriate here. We have modified\nthis sentence in the revised manuscript. We have also added more descriptions of AUC to the\nrevised manuscript for better understanding (lines 193-194, 198-199, 388-395).\n11. Line 182: The phrase “four times the training time” is misleading without context. This\nscaling is due to the use of an ensemble of size 4, which does not appear to be\nmentioned earlier. It would be more accurate to state that training time scales linearly\nwith the number of models in the ensemble. Similarly, inference time for both\nensembles and MC dropout (a pseudo-ensemble) also scales linearly with the\nensemble size—here, 4.\nResponse:\nThank you. We have revised the manuscript accordingly (lines 208-218).\n12. Line 186: The authors mention that expectation-maximization (EM) increases\ninference cost for the GMM model. However, EM is typically used during GMM\ntraining, not inference. At prediction time, GMM inference should be relatively\nlightweight. Could the authors clarify what is contributing to the increased cost of\nGMM inference, and why it is reported as slower than an ensemble?\nResponse:\nThank you for pointing this out. In the original figure, we mistakenly included the EM time as\npart of the overall inference time, while it should actually be attributed to the training process.\nIn fact, GMM inference is indeed computationally lightweight. In the revised figure, we have\nused a light blue dashed bar to represent the EM extra time and corrected the inference time\nfor GMM accordingly.\n13. Figure 3c: Could the authors consider reporting the force RMSE instead of the MAE?\nFor vector quantities, the L1 norm is not invariant under rotation, so MAE can yield\ninconsistent values depending on molecular orientation. The RMSE (L2 norm) does\nnot suffer from this problem and would be more appropriate here.\nResponse:\nThank you for your suggestion. We have computed the RMSEs, and both MAE and RMSE\nare presented in the revised figure. As can be seen, the results of RMSE are consistent with\nMAE.\n14. Figure 4c: Up to this point, the authors emphasize that their model yields directional\nuncertainties. However, in Figure 4c, the y-axis is labeled simply as uncertainty. Since\nthe model appears to predict a tensorial uncertainty per configuration, it would be\nhelpful to explain how this scalar value is derived from the full prediction. A brief\nexplanation or a pointer to the relevant equation in the Methods section (where it is\nexplained) would clarify this.\nResponse:\nThank you. We have added a reference to Equation 7 (the uncertainty calculation formula) to\nthe caption of Figure 4.\n15. Lines 200–203: This brief explanation of how directional uncertainty helps identify\nchemical reactivity is insightful. Consider moving a version of this explanation to the\nIntroduction, where it could better motivate the relevance of the method.\nResponse:\nThank you. We have revised the manuscript accordingly (lines 101-102).\n16. Lines 206–207: The phrasing “standard MD and eIP-UDD” is potentially confusing.\nSince both simulation types are carried out using the trained eIP model, I suggest\nrephrasing to something like “we performed both standard MD and UDD simulations\nusing our eIP model”. Alternatively, if the term eIP-UDD is retained throughout the\nmanuscript, it should be clearly defined at first occurrence to avoid the impression\nthat two different force fields are being compared.\nResponse:\nThank you. We have clearly defined the term eIP-UDD and revised the manuscript\naccordingly (lines 238-242).\n17. Line 215: The statement that the simulations are “stable after 50 ps” might be\nmisleading. Do the authors mean that the simulations remain stable for the full 50 ps?\nAlso, the current phrasing does not communicate what advantage the UDD\nsimulation provides, especially since both UDD and standard MD lead to stable\nmodels after two active learning iterations. Could the authors clarify the added value\nof UDD in this context?\nResponse:\nWe apologize for the confusion caused by the original expression. Yes, we mean that the\nsimulations remain stable for the full 50 ps. This experiment was intended as a deliberate,\nmulti-stage demonstration of the eIP-UDD simulation's capabilities. Our goal was to\nshowcase how UDD behaves under different sampling conditions: the initial simulation\nrevealed the collapse due to insufficient sampling; the second highlighted UDD's\neffectiveness in exploring and acquiring higher-uncertainty configurations to improve stability;\nand the third demonstrated its behavior after the local configurational space had been\nsufficiently sampled.\nWe did not adjust UDD parameters throughout all iterations for comparability. The similar\nbehavior of UDD and MD simulations indicates that the corresponding local configurational\nspace has been sufficiently sampled under the given parameters. We would like to emphasize\nthat these distinct outcomes at each stage are fully tunable. As we demonstrated in the\nrevised Supplementary Section S7 and Figures S9-S12, a broader configurational space can\nbe further explored by adjusting the UDD parameters. We have revised the analysis of Figure\n4 in the revised main text (lines 245-254) to clarify this point.\n18. Figure 5: The caption should specify which data set is used. In panels (d) and (g), the\nabrupt changes in potential energy at the start of the trajectories are concerning.\nHowever, they could plausibly be attributed to the lack of proper thermalization. Was\nthere a heating or equilibration stage prior to the MD runs shown here? If not, a brief\nexplanation of why this was omitted—or a clarification in the figure caption—would\nhelp interpret the results.\nResponse:\nThere was no heating or equilibration stage prior to the MD runs shown here. In panel (g), the\ninitial configuration of PDMS was constructed in-house and not pre-relaxed. The decrease in\npotential energy corresponds to the structural relaxation process. In panel (d), the initial\nconfiguration of LiFePO4 was sourced from the Materials Project and had already been\noptimized using DFT. Therefore, no significant decrease in potential energy is observed at the\nbeginning of the simulation. The brief rise in potential energy during the early stage of the\nUDD simulation occurs due to modifications made to the potential energy surface, which allow\nthe system to escape the original local minimum. The discussion regarding potential energy\nchanges is added to the revised manuscript (lines 268-275). The dataset used for training the\nuniversal potential is included in the revised caption.\n19. Lines 228 and 333–334 (configurational entropy): A reference to the relevant Methods\nor SI section would be helpful here. Also, is there a standard citation for the specific\napproach used? Furthermore, could the authors expand on how configurational\nentropy is computed in this work? For example, is a long unbiased MD simulation first\nperformed to collect histograms of key order parameters, which are then used to\nestimate configuration probabilities? Or are the probabilities inferred from the training\nset? If so, this method relies on the assumption of a sufficiently unbiased and ergodic\nsampling, which may not hold in the case of data sets obtained through active\nlearning. In such scenarios, the resulting histograms may not reflect a proper\nBoltzmann distribution, making entropy estimates difficult to interpret.\nResponse:\nConfigurational entropy is an important statistical measure in molecular dynamics simulation\nstudies. The methods of estimating configuration entropy are reported in literature [1,2]. Here,\nwe adopted an approximate calculation method. The configurational entropy is defined as the\nShannon entropy [3,4]:\n𝑆 =−∑𝑝 ln𝑝,\n𝑖 𝑖\n𝑖\nwhere 𝑝 denotes the probability of the system being in state i. Next, we projected states\n𝑖\nonto a discretized order parameter grid and then calculated the frequency of these order\nparameters within a simulation trajectory. Strictly speaking, this does not represent true\nconfigurational entropy, as the trajectory length may fall short of meeting the requirement of\nergodicity. The purpose of this calculation was to evaluate the configurational diversity within\nthe simulation trajectory.\nThe detailed procedure for configurational entropy calculation is elaborated and explicitly\nreferenced in the revised manuscript (lines 265-267, 398-406).\n[1] Karplus, M. & Kushick, J. N. Method for estimating the configurational entropy of\nmacromolecules. Macromolecules 14, 325–332 (1981)\n[2] Peter, C., Oostenbrink, C., Van Dorp, A. & Van Gunsteren, W. F. Estimating entropies from\nmolecular dynamics simulations. The Journal of Chemical Physics 120, 2652–2661 (2004)\n[3] Baxa, M. C., Haddadian, E. J., Jha, A. K., Freed, K. F., & Sosnick, T. R. Context and Force\nField Dependence of the Loss of Protein Backbone Entropy upon Folding Using Realistic\nDenatured and Native State Ensembles, Journal of the American Chemical Society 134,\n15929–15936 (2012)\n[4] Klyshko, E. et al. Functional protein dynamics in a crystal. Nature Communications 15,\n3244 (2024)\n20. Line 258: When referring to \"residue\", do the authors mean \"residual\" as in the\ndifference between the predicted and reference values?\nResponse:\nYes, it should be \"residual\". Thank you for pointing out the error. It has been corrected.\n21. Equation 7: Have the authors considered using alternative aggregation functions for\nestimating the global uncertainty, beyond the mean? For instance, in larger systems\nlike silica, it is possible that a few atoms experience large force prediction errors,\nwhich could lead to failure in the simulation—even though the average uncertainty\nremains low. A maximum or weighted quantile might better reflect risk in such edge\ncases. Also, should both sides of the equation include a square term? As written, it\nappears the square of the per-atom uncertainties may be missing.\nResponse:\nWe acknowledge that taking the mean may result in local information loss. However, selecting\nthe maximum value may also be problematic due to the intrinsic errors in uncertainty\nquantification. The weighted quantile should be a better choice, but the specific weights need\nto be adjusted. This has been discussed in the revised manuscript (lines 322-325).\nIn addition, there was a typo in Equation 6: each term in the equation should not have a\nsquare. Accordingly, Equation 7 also does not include squares. This is consistent with\nEquations S9 and S10.\n22. Equations 8 and 9: The motivation and derivation for these equations could be\nexpanded. As it stands, they are difficult to follow for readers unfamiliar with Bayesian\nquantile regression. A brief explanation of how these formulations emerge from\nthe underlying assumptions or loss functions would improve accessibility.\nResponse:\nWe have expanded the derivations for these equations (lines 326-332) to make them more\naccessible and easier to follow.\n23. Lines 250 / 351: Is it advisable to start simulations directly at 300 K, particularly in\nearly active learning iterations when the MLIP may still be unreliable? This practice\ncould increase the likelihood of simulation instabilities or system blow-ups. Could the\nauthors comment on this choice and whether a lower initial temperature or a\ntemperature ramp was considered?\nResponse:\nBased on our experience, we recommend starting simulations directly at 300 K. We have\nobserved that when the training data is insufficient, even if the temperature is reduced, the\nsimulation tends to collapse after a certain period. Although the duration before collapse may\nbe extended, the dynamics also slow down at lower temperatures, which may not be very\nhelpful in sampling new data.\nIn practice, we are usually able to construct a reasonably good dataset in the initial stage. For\nexample, we can run multiple simulations at various temperatures and pressures using\nempirical force fields or universal potentials, and then sample from each trajectory and\nperform DFT calculations. At this stage, simulation stability is generally not a major issue, but\nthe simulation results may not be correct. Therefore, multiple rounds of active learning\niterations are needed to refine the training set.\n24. Very minor recommendation regarding SI front page: It would be helpful to include a\nshort title page for the SI, stating the title of the paper, author names, and affiliations.\nThis is standard in many journals and aids in readability and citation.\nResponse:\nThank you. A title page for SI has been added.\n25. Section S3: The sentence “When ablating locality, we chose global embedding to\nestimate the uncertainty associated with the energy.” is unclear. There is no “global\nembedding” in the PaiNN architecture. If the authors refer to pooling features across\natoms or to removing the spatial locality constraint in uncertainty prediction, they\nshould reformulate this more precisely.\nResponse:\nYes, we refer to pooling features across atoms. It has been corrected accordingly.\n26. Section S4: “In contrast to the standard PaiNN model, the eIP model incorporates\nadditional fully connected layers.” It is unclear whether the authors refer to changes in\nthe output head or in the message passing layers themselves. The original PaiNN\nalready defines multiple output architectures, such as simple feed-forward networks\nfor scalar outputs and gated equivariant blocks for vector-valued quantities. If the\nauthors modified or replaced these, it would be valuable to include a visual diagram\nof the architecture, especially to clarify whether these changes are post-message-\npassing or occur in parallel.\nResponse:\nWe have added a figure illustrating the architecture of the eIP output layer (Figure S7).\n27. SI Eq. 5 and Eq. 7: Could the authors specify the dimensionality or size of the term ε?\nResponse:\nIt is a scalar and will be broadcast to an Nx3 shape during addition according to Python's\nbroadcasting rules. This has been clarified accordingly.\n28. SI Eq. 9 vs Main Text Eq. 6: There appears to be an inconsistency between Eq. 9 in\nthe SI and Eq. 6 in the main text. In the SI, the variances are summed directly, while\nin the main text, the squares of the variances are added before taking the square\nroot. The authors should clarify this discrepancy and ensure consistency between\nformulations.\nResponse:\nWe apologize for the typo in the main text equation 6: each term in the equation should not\nhave a square. This is consistent with Equation S9.\n29. SI, page 9 (\"final layer\"): The term “final layer” could refer to the last layer of the\nmessage passing network or to the last layer of the output head. Please clarify which\nis meant here, especially as this distinction may affect how uncertainty is computed\nor aggregated.\nResponse:\nWe mean the last layer of the message passing network. It has been clarified accordingly.\nThank you.\nReviewer 2\nThis report details the use of evidential deep learning to capture epistemic uncertainty in\nMLIPs at reduced computational cost as compared to standard methods such as ensembling,\nMonte Carlo dropout, GMM, and MVE. Uncertainty-driven dynamics simulations are used\nalongside their eIP method in an active learning process. Their eIP method is of significance\nto the field of MLIPs and will help to model epistemic uncertainties without the need for\nexpensive model ensembling.\n1. In the discussion of aleatoric versus epistemic uncertainty, while it is true that strict\ncuration of a dataset would limit the aleatoric uncertainty, this type of careful data\ncuration is often not what is done in practice, especially in newer foundation\nMLIPs that often combine multiple large datasets. It should be noted in the main text that\naleatoric uncertainty is present in many pretrained foundation/universal models. In\nparticular, this should be discussed for the MPtrj dataset, as there are known sources\nof aleatoric uncertainty in that dataset. See this report for example:\nhttps://www.nature.com/articles/s41524-025-01572-y\nResponse:\nThank you for the insightful comment and the helpful reference. We agree that many MLIPs\nincorporate multiple datasets of varying accuracy for training. In practice, obtaining high-\naccuracy labels through rigorous ab initio calculations is highly challenging. As mentioned in\nthe paper you provided, aleatoric uncertainty may arise from inconsistent Hubbard U\ncorrection, varying convergence criteria, and so on.\nWe have revised the manuscript (lines 120-124) to discuss this. We have also added\nadditional experiments to the revised SI (S2 Aleatoric uncertainty) to provide a more\ncomprehensive discussion on the origins and implications of aleatoric uncertainty for MLIPs.\n2. Model ensembling is known to produce overconfident uncertainties during inference.\nThe authors use an evidence regularization function to reduce model overconfidence.\nMore discussion should be given concerning this topic. A demonstration should be\nincluded to show the effect of the evidence regularization function on reducing\noverconfident uncertainties.\nResponse:\nThis evidence regularization function was initially proposed by Amini et al. [1],\n𝐿𝑅 =|𝑦 −𝔼[𝜇 ]|⋅Φ=|𝑦 −𝛾|⋅(2𝜈+𝛼),\n𝑖 𝑖 𝑖 𝑖\nwhere the total evidence is defined as Φ=2ν+α . In our work, we used a modified version\nintroduced by Huttel et al. [2], which replaced the MAE loss with a tiled loss, and the total\nevidence is redefined as Φ=2ν+α+1/β.\nThe total evidence Φ is inversely correlated with predictive uncertainty. When predictions are\ninaccurate, the model learns to reduce its confidence (i.e., output lower values for 𝜈 and 𝛼,\nor a higher value for 𝛽). Consequently, this regularization term effectively mitigates\noverconfidence.\nTo evaluate the efficacy of the regularization term, we conducted an ablation experiment by\nsetting the regularization weight 𝜆 to 0, 0.01, 0.1, and 10, respectively. We then quantified\nthe degree of overconfidence using the miscalibration area below the 𝑥 =𝑦 line [3]. We also\ncompared it with the ensemble method. The results are shown in Table S1. We have revised\nthe manuscript (lines 333-339) and added a discussion to the Supplementary section S4.\n[1] Amini, A., Schwarting, W., Soleimany, A. & Rus, D. Deep evidential regression. Advances\nin Neural Information Processing Systems 33, 14927–14937 (2020).\n[2] Huttel, F. B., Rodrigues, F., & Pereira, F. C. Deep evidential learning for bayesian quantile\nregression. arXiv preprint arXiv:2308.10650 (2023).Bayesian\n[3] Vita, J. A., Samanta, A., Zhou, F. & Lordi, V. Ltau-ff: loss trajectory analysis for uncertainty\nin atomistic force fields. Machine Learning: Science and Technology 6, 015048 (2025)\n3. Because the silica glass dataset used in this work is not commonly used in\nbenchmarking MLIPs, more information is needed describing the dataset, including\nthe differences between the ID and OOD datasets partitioned in reference 30.\nResponse:\nThank you for your suggestion. We have added a detailed description of this dataset in the\nrevised manuscript (lines 350-360).\n4. Please provide the empirical force field used to generate the classical MD water\ndataset used during eIP-UDD.\nResponse:\nThank you for the reminder. We used the SPC/E force field. This has been clarified in the\nrevised manuscript (line 362).\nReviewer 3\nIn the manuscript “Evidential Deep Learning for Interatomic Potentials” the authors Han Xu et\nal describe a method for training machine learning interatomic potentials to predict\nparameters in a distribution rather than a singular energy and force value. This particular\nparameters inferenced in this distribution are then used to indicate whether a specific atomic\nconfiguration in active learning should be simulated with electronic structure theory or not.\nThe authors due a commendable job of comparing this method to other uncertainty\nquantification methods and link it to previously developed accelerated sampling techniques\n(UDD). In general, I find the manuscript to be well can clearly written and provide a valuable\nmethodology for reducing the computational cost associated with active learning. My only\nconcern is that the technique demonstrated here is only of use to those in the community\nactively building MLIP training data with UQ methods. This is a relatively small community\nand thus this paper may not be suitable for nature communications simply because of the\nniche it occupies. This is a decision to be made by the editor.\n1. My only real objection in the manuscript is that the ensemble method is shown to be a\nclear factor of 4 more expensive for training and inference tasks. This is not quite\ncorrect for various practical and technical reasons. For training, most HPC nodes\nhave multiple GPUs. Due to the limitation of data transmission/memory\nsynchronization, training 4 separate models on 4 separate GPUs only takes about 2\ntimes longer than training 1 model on 4 GPUs. This is obviously package/model\ndependent but would guess similar trends would hold for other Pytorch based\nmodels. Further, in inference, ensemble-based methods, when properly implemented,\nare not “# of networks” more expensive than a single model. This is because things\nlike mapping of atoms to symmetry functions (which are shared between models) can\nbe re-used for each model. The relatively cheap tensor operations of the MLIP are the\nonly things that need to be repeated for each element of the ensemble. This\npractically saves 40% of calculation time in the ensemble, making it much more\ncompetitive than Fig 3 suggests. I realize these details are difficult to discuss in a\ngeneral sense, but possibly a footnote on ensemble timings is worth\nincluding.\nResponse:\nWe greatly appreciate your insightful comment. We acknowledge that our statement about the\nensemble method being a factor of 4 more expensive is an oversimplification that does not\naccount for technical optimization. We agree that a more detailed discussion is necessary to\naccurately describe the computational cost. We have revised the manuscript (lines 208-216)\nto clarify this point.\n2. This manuscript also suggests interesting questions outside of the present scope:\nDoes aleatoric uncertainty actually relate to noise in the data labels? This is easy to\nexplore in an electronic structure context because the convergence threshold of\nbasically every calculation is adjustable. It would be interesting to numerically verify\nthat the author’s label of aleatoric uncertainty goes up when lower convergence\nthresholds are used. Further, does the ability to train a specific aleatoric uncertainty\nparameter enable the model to more efficiently use poorly converged simulations,\npotentially saving time on electronic structure calculations?\nResponse:\nThis question is very interesting and important. In response, we conducted several additional\nexperiments. We agree that a numerical experiment to study the relationship between the\nconvergence threshold and aleatoric uncertainty would be ideal. However, performing such an\nanalysis on the silica glass dataset, which contains configurations with 699 atoms, would\nrequire several months of computational time. Furthermore, label noise can also be\nintroduced by inconsistencies in the choice of basis set, functional, and Hubbard U correction\n[1]. As a practical alternative, we used the universal potential trained in this work on the MPtrj\ndataset as a proxy for low-accuracy labels. The difference between the original DFT forces\nand the low-accuracy forces is shown in Figure S4.\nSurprisingly, the model trained on the low-accuracy labels exhibits significantly lower\nepistemic uncertainty than the model trained on the original DFT-labels, while both models\nhave similar aleatoric uncertainty. To understand this result, we hypothesize that the\nuncertainty values obtained from different models may not be directly comparable, as these\nvalues can be influenced by the regularization term in the loss function. However, the\nmagnitudes of the two types of uncertainties derived from the same model should be\ncomparable and yield more reasonable conclusions. For the low-accuracy model, the mean\naleatoric uncertainty (0.455) and the mean epistemic uncertainty (0.542) are similar. In\ncontrast, for the high-accuracy model, the mean aleatoric uncertainty (0.482) is only 1/200 of\nthe mean epistemic uncertainty (96.519). The overwhelming dominance of epistemic\nuncertainty in the high-accuracy model indicates that its total uncertainty is primarily reducible\nwith more data, whereas the uncertainty in the low-accuracy model is constrained by\nirreducible noise.\nRegarding practical application, aleatoric uncertainty can be useful for active learning.\nAlthough we currently lack a clear strategy for efficiently using poorly converged data, it is\npossible to filter out data with high aleatoric uncertainty but low epistemic uncertainty to\nimprove the consistency of labels. This could be particularly advantageous for training\nfoundation MLIPs that combine multiple large datasets.\nWe acknowledge that due to computational constraints, our experimental validation of this\napproach remains limited. A thorough investigation into the practical application of aleatoric\nuncertainty is a critical direction for our future work.\nThe additional experimental results and discussions have been added to the revised SI,\nsection S2 \"Aleatoric uncertainty\".\n[1] Bilbrey, J.A., Firoz, J.S., Lee, MS. et al. Uncertainty quantification for neural network\npotential foundation models. npj Comput Mater 11, 109 (2025)"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I thank the authors for their diligence in addressing all my concerns and performing the additional experiments. I recommend\nthe article for publication in its current form.", "reviewer_2": "I appreciate the authors' careful consideration of the reviewers' comments and believe the article is ready for publication.\nA knowledgeable researcher in the field could run the code, but a novice may have trouble. A more detailed README file\nwould help.", "reviewer_3": "The authors have done a thorough job responding to reviewer feedback. I have no objections to this manuscript."}, "author_response": ""}]

Decision letter