GPR_calculator: An On-the-Fly Surrogate Model to Accelerate Massive Nudged Elastic Band Calculations

We present GPR_calculator, a package based on Python and C++ programming languages to build an on-the-fly surrogate model using Gaussian Process Regression (GPR) to approximate expensive electronic structure calculations. The key idea is to dynamically train a GPR model during the simulation that can accurately predict energies and forces with uncertainty quantification. When the uncertainty is high, the expensive electronic structure calculation is performed to obtain the ground truth data, which is then used to update the GPR model. To illustrate the power of GPR_calculator, we demonstrate its application in Nudged Elastic Band (NEB) simulations of surface diffusion and reactions, achieving 3-10 times acceleration compared to pure ab initio calculations. The source code is available at https://github.com/MaterSim/GPR_calculator.

💡 Research Summary

The manuscript introduces GPR_calculator, a hybrid Python‑C++ package that integrates Gaussian Process Regression (GPR) as an on‑the‑fly surrogate for expensive electronic‑structure evaluations within Nudged Elastic Band (NEB) calculations. The authors recognize that density functional theory (DFT), while accurate, becomes a computational bottleneck when hundreds to thousands of energy and force evaluations are required to converge a minimum‑energy path (MEP). To alleviate this, GPR_calculator couples a conventional “base” calculator (any DFT code interfaced through ASE) with a GPR model that predicts energies and forces for each NEB image and simultaneously provides a variance‑based uncertainty estimate. When the predicted uncertainty exceeds a user‑defined threshold, the package automatically invokes the base calculator, obtains the true DFT data, and updates the GPR model in real time. This adaptive scheme ensures that only the most informative configurations are evaluated with the expensive method, dramatically reducing the total number of DFT calls.

The technical core rests on a rotation‑invariant SO(3) descriptor derived from the power spectrum of spherical‑harmonic expansions of the local atomic density. Each atom is represented by a 30‑dimensional vector (n_max = 3, l_max = 4) that captures both radial and angular information within a cutoff of 5–6 Å. The descriptor’s analytical derivatives with respect to atomic coordinates are also computed, enabling the construction of force kernels. Two kernel families are offered: a radial‑basis‑function (RBF) kernel and a dot‑product kernel, both built on a cosine similarity metric d(x₁,x₂) that can be raised to an integer power ζ to sharpen discrimination. The kernel matrix C is formed from energy‑energy (k_ee), energy‑force (k_ef), and force‑force (k_ff) blocks, and the GPR hyper‑parameters (signal variance θ, length scale l, power ζ, and noise variance β) are optimized by maximizing the log‑likelihood using the L‑BFGS‑B algorithm.

Computationally, the authors mitigate the O(N³) cost of inverting C by employing a low‑rank approximation (truncating small eigenvalues) and Cholesky factorization, reducing the scaling to O(N M²) where M ≪ N. When a new training point is added, the covariance matrix is expanded via block augmentation, and the existing factorization is updated efficiently. Parallelization is achieved at the image level and within kernel evaluations using MPI/OpenMP, allowing the package to exploit modern multi‑core clusters.

Performance is demonstrated on two prototypical surface processes: adatom diffusion on Pt(111) and a CO desorption reaction on a catalytic surface. Compared with pure DFT‑driven NEB, GPR_calculator achieves a 3–10× speed‑up (average 4.2×) while maintaining saddle‑point energy errors below 0.02 eV. The authors benchmark against earlier machine‑learning approaches, notably neural‑network MLFFs and static GPR models, highlighting that the on‑the‑fly uncertainty‑driven sampling dramatically improves data efficiency and eliminates the need for a separate validation step after training.

Limitations are acknowledged: the GPR framework still scales poorly beyond a few thousand training points, suggesting future integration of sparse GPR, local kernels, or deep kernel learning. The descriptor calculation, though compact, adds overhead for very large systems, motivating exploration of lighter alternatives such as SOAP or graph‑based embeddings. Nonetheless, the authors argue that the methodology is readily extensible to other transition‑state theories, metadynamics, and free‑energy surface constructions.

In conclusion, GPR_calculator provides a practical, open‑source solution for accelerating NEB simulations by marrying rigorous uncertainty quantification with adaptive surrogate modeling. Its ability to automatically balance accuracy and computational cost positions it as a valuable tool for high‑throughput catalyst design, materials discovery, and any workflow where transition‑state searches dominate the computational budget.