Source-linked AI summary
Active learning of reactive Bayesian force fields: Application to heterogeneous hydrogen-platinum catalysis dynamics
Jonathan Vandermause, Yu Xie, Jin Soo Lim, Cameron J. Owen, Boris Kozinsky
TL;DR
Reactive force-field modeling must balance the accuracy needed for chemical reactions with the computational cost and expertise required by existing approaches. FLARE++ trains a reactive many-body force field on the fly using Bayesian uncertainty and then maps it to a faster polynomial model. For hydrogen reactions on Pt(111), the resulting model is twice as fast as recent Pt/H ReaxFF while being considerably more accurate, and the training requires three days of wall time.
Problem
Accurate reactive force fields are difficult to obtain because ab initio methods are expensive and flexible alternatives such as ReaxFF require substantial time, effort, and expertise to parameterize.
Method
FLARE++ uses Bayesian uncertainty from a sparse Gaussian process during molecular dynamics to select ab initio data, then maps the trained model to an equivalent polynomial model with training-set-size-independent prediction cost.
Results
On Pt(111) hydrogen splitting and recombination, the trained model is twice as fast as a recent Pt/H ReaxFF force field and considerably more accurate.
Takeaways & Limitations
The method reduces the time and effort needed to train fast and accurate reactive force fields for complex systems.
Abstract
from arXiv · showhide
Accurate modeling of chemically reactive systems has traditionally relied on either expensive ab initio approaches or flexible bond-order force fields such as ReaxFF that require considerable time, effort, and expertise to parameterize. Here, we introduce FLARE++, a Bayesian active learning method for training reactive many-body force fields on the fly during molecular dynamics (MD) simulations. During the automated training loop, the predictive uncertainties of a sparse Gaussian process (SGP) force field are evaluated at each timestep of an MD simulation to determine whether additional ab initio data are needed. Once trained, the SGP is mapped onto an equivalent and much faster model that is polynomial in the local environment descriptors and whose prediction cost is independent of the training set size. We apply our method to a canonical reactive system in the field of heterogeneous catalysis, hydrogen splitting and recombination on a platinum (111) surface, obtaining a trained model within three days of wall time that is twice as fast as a recent Pt/H ReaxFF force field and considerably more accurate. Our method is fully open source and is expected to reduce the time and effort required to train fast and accurate reactive force fields for complex systems.
I. INTRODUCTION
Reactive molecular dynamics needs force fields that can describe bond making and breaking while remaining accurate and inexpensive enough for long simulations. FLARE++ addresses the training challenge with on-the-fly Bayesian active learning and an accelerated polynomial model, demonstrated for hydrogen reactions on Pt(111).
- Motivation: Reactive MD requires a potential energy surface that is accurate for chemical bond breaking and formation yet cheap enough for long-timescale simulations.Ab initio methods are difficult to extend beyond a few hundred atoms or a few hundred picoseconds, while fixed-bond force fields cannot describe reactions.
- Motivation: Manual machine-learning force-field training is especially difficult for reactive systems because relevant transition paths and sampling requirements are unknown in advance.Such training can require months of effort, substantial expertise, and significant computing resources.
- FLARE++ method: FLARE++ uses Bayesian uncertainties from a sparse Gaussian process force field during molecular dynamics to select structures requiring additional ab initio calculations.The model is updated when uncertainty exceeds a chosen threshold, enabling automated on-the-fly training.
- FLARE++ method: The trained sparse Gaussian process is mapped to an equivalent polynomial model whose prediction cost is independent of training-set size.This mapping targets the linear sparse-set cost of kernel-based predictions while preserving the model's learned behavior.
- Application: For hydrogen splitting and recombination on Pt(111), the accelerated model is more than twice as fast as a recent ReaxFF force field and considerably more accurate.The model gives an activation energy for hydrogen turnover in close agreement with the reference.
II. RESULTS
FLARE++ combines Bayesian uncertainty-guided active learning with a polynomial mapping that preserves many-body force-field behavior while making prediction cost independent of sparse-set size.
- Active learning: The method automatically constructs training and sparse sets for reactive many-body force fields during molecular dynamics.Bayesian uncertainties determine when DFT data are required, while novel environments are selectively added to the sparse set.
- Active learning: When local-energy uncertainty exceeds the DFT tolerance, the simulation pauses for a DFT calculation whose energy, forces, and stresses update the model.Predictions below tolerance are accepted for the next MD step.
- Many-body descriptors: The descriptor maps local environments to symmetry-preserving many-body features that encode angular information through three-body contributions.The resulting invariant descriptor is used as input to the sparse Gaussian process.
- Acceleration: The trained sparse Gaussian process is mapped to an equivalent polynomial model whose prediction cost is independent of sparse-set size.An integer kernel power determines polynomial order and the learned force field’s body order.
- Model selection: For Pt/H, ξ = 2 has much higher likelihood than ξ = 1 and nearly the same likelihood as ξ = 3, so the authors select a five-body model.Likelihood decreases for ξ > 3.
- Acceleration: The quadratic implementation achieves nearly twenty-fold acceleration over standard sparse-GP prediction and exceeds the speed of recent Pt/H ReaxFF by more than twofold.The accelerated model enables more efficient machine-learning reactive MD than the classical comparison force field.
B. On-the-fly training of a reactive Pt/H force field
On-the-fly training learns a reactive Pt/H force field from uncertainty-triggered DFT calculations, then validates it on interpolation, extrapolation, and materials properties.
- Training procedure: The Pt/H training simulation used 216 DFT calls over 3.7 ps and produced nearly 50,000 energy, force, and stress labels.Most DFT calls occurred during the reactive Pt/H simulation rather than the single-element runs.
- Training procedure: Energy predictions agreed with DFT to within 1 meV/atom during the Pt/H training simulation.The simulation began with gas-phase H2 and a hydrogen-covered Pt surface at 1500 K to facilitate rare recombination events.
- Active learning: A hydrogen-bond formation event triggered two DFT calls, demonstrating that the active-learning loop detects novel environments and adds them to training.The event occurred during the first recombination sequence at approximately 0.3 ps.
- Validation: The final model was assembled from four independent simulations and evaluated on properties both represented and not explicitly included in training.Validation covered energies, forces, stresses, bulk platinum, hydrogen, and Pt/H adsorption behavior.
- Validation: The model predicts the fcc-Pt lattice constant within 0.1% and bulk modulus and elastic constants within 6% of DFT, while ReaxFF overestimates C44 by nearly 200%.This comparison concerns bulk-platinum properties.
C. Large-scale reactive MD
The accelerated reactive force field supports large-scale Pt(111) simulations of hydrogen splitting and recombination across temperatures, yielding activation energies consistent with experiment.
- Simulation setup: The trained SGP was mapped to an accelerated quadratic model for fixed-temperature simulations of hydrogen splitting and recombination on Pt(111).The simulations used a six-layer 12-by-12 slab, 864 Pt atoms, and 448 H atoms.
- Reaction dynamics: After equilibration, hydrogen splitting and recombination rates became roughly equal, while higher temperatures produced lower equilibrium surface coverage.Initially, recombination exceeded adsorption as the surface approached equilibrium coverage.
- Reaction kinetics: The estimated activation energy was 0.25(2) eV, and a doubled-vacuum test gave 0.20(3) eV.Both estimates were obtained from reaction-rate fits over the final 300 ps of each simulation.
- Reaction kinetics: Both simulated activation-energy estimates agree with the experimentally measured value of 0.23 eV.The agreement was observed across the standard and doubled-vacuum simulations.
III. DISCUSSION
The authors present FLARE++ as a unified route to accurate, efficient reactive force fields that reduces training burden and broadens reactive MD applications.
- Discussion: FLARE++ achieves excellent accuracy relative to DFT while remaining computationally competitive with ReaxFF for reactive molecular dynamics.The discussion frames this as a combined accuracy and efficiency result for reactive force fields.
- Discussion: The method is intended to reduce the time, effort, and expertise needed to train accurate reactive force fields, producing a high-quality model within days with minimal supervision.This goal addresses the practical burden of parameterizing reactive force fields.
- Discussion: FLARE++ bridges GAP-like Bayesian uncertainty with the computational efficiency of parametric force fields such as SNAP, qSNAP, and MTP.The approach combines uncertainty information from kernel-based models with efficient parametric prediction.
- Outlook: The authors expect the unified approach to extend reactive MD toward complex systems that have remained difficult to model computationally.They identify biochemical reactions and more complex heterogeneous catalysts as areas of interest.
IV. METHODS
The methods implement sparse Gaussian-process force fields, accelerated polynomial mappings, and molecular-dynamics, machine-learning, and DFT calculations.
- The methods cover SGP force fields, polynomial acceleration, and MD, ML, and DFT computational details.These components are presented as the paper’s key methodological elements.
- Three-dimensional arrows denote spatial vectors, while bold symbols denote vectors and tensors in general feature spaces.
A. Sparse Gaussian process force fields
The force field expresses total energy through atom-centered local environments, symmetry-aware descriptors, and a kernel measuring environmental similarity.
- Local environments are mapped to descriptor vectors, and a kernel quantifies similarity between pairs of environments.This provides the model inputs and similarity measure used by the force field.
- The SGP represents total potential energy as a sum of local energies assigned to atom-centered environments.Each environment contains neighboring species and interatomic vectors within species-dependent cutoff spheres.
- Species-dependent cutoffs provide additional flexibility, with shorter H-H and H-Pt cutoffs than the Pt-Pt cutoff in the Pt/H system.The passage attributes this cutoff choice to the requirements of the Pt/H system.
2. Describing local environments
Local environments are encoded with symmetry-preserving many-body descriptors built from radial functions and spherical harmonics, then reduced to unique invariant components.
- The descriptor construction produces many-body features satisfying rotational, permutational, translational, and mirror symmetry.Rotational invariance is obtained by converting covariant descriptors into invariant ones.
- Interatomic vectors are expanded with radial basis functions and real spherical harmonics, using a cutoff that smoothly vanishes at the cutoff radius.The radial basis is defined on [0, 1], and the cutoff depends on scaled interatomic distance.
- A covariant tensor is formed by summing basis functions over neighboring atoms by species, then contracted to obtain a rotationally invariant descriptor.The invariant descriptor uses the spherical-harmonic sum rule and removes redundancy from interchangeable indices.
3. Making model predictions
Predictions combine sparse-environment kernels with Bayesian uncertainty estimates for local, total, surface, and binding energies.
- The SGP predicts local energy as a weighted sum of kernels between a query environment and representative sparse environments.The kernel is a normalized dot product raised to integer power ξ, analogous to the SOAP kernel.
- The model uses training labels for potential energies, forces, and virial stresses, with noise values represented in a diagonal matrix.The final model uses force, energy, and stress noise hyperparameters, and QR decomposition avoids unstable direct inversion.
- Predictive variances extend from local energies to total, surface, and binding energies through covariance combinations under the DTC approximation.The resulting confidence regions are used for the surface and binding energies reported in Fig. 3.
- A scaled local-energy variance lies between 0 and 1 and is independent of kernel hyperparameters, providing the uncertainty measure used to guide active learning.This measure is used to identify uncertainty in local environments during the training protocol.
5. Optimizing hyperparameters
The SGP hyperparameters are optimized by maximizing the DTC log marginal likelihood, balancing model complexity against fit quality during on-the-fly training.
- The DTC log marginal likelihood balances model complexity against the quality of fit when selecting hyperparameters.Its first term penalizes complexity, while its second measures fit quality.
- After the first Nhyp SGP updates, kernel hyperparameters σ, σE, σF, and σS are optimized with L-BFGS using gradients of L.
- The log marginal likelihood is also used to evaluate different hyperparameter choices during on-the-fly runs.
B. Mapping to an equivalent polynomial model
The trained SGP can be rewritten as an equivalent polynomial model in local-environment descriptors, eliminating loops over sparse environments and making prediction cost independent of training-set size.
- Once the sparse-set terms are gathered into a rank-ξ tensor β, SGP mean predictions avoid loops over sparse points and accelerate substantially.This reformulation is especially useful for small ξ.
- For ξ = 1, mean predictions become linear in the descriptor and require a single dot product for local-energy evaluation.
- The polynomial model’s prediction cost is independent of the number of sparse environments, unlike standard SGP prediction.
- For ξ = 2, mean predictions are quadratic in the descriptor and evaluate through a vector-matrix-vector product.
- The implementation uses ACE descriptors and provides code for training SGPs and mapping them onto accelerated quadratic models.