PaperlayerView exact revision on arXiv

Source-linked AI summary

On-the-fly machine learning force field generation: Application to melting points

Ryosuke Jinnouchi, Ferenc Karsai, Georg Kresse

arXiv:1904.12961v3cond-mat.mtrl-sciphysics.comp-ph

TL;DR

Computationally tractable and broadly applicable first-principles melting-point prediction remains difficult, while conventional machine-learning force fields demand large datasets and extensive supervision. This paper develops an on-the-fly Bayesian force-field generation method integrated into an electronic-structure code, triggering first-principles calculations only for uncertain configurations. Applied to Al, Si, Ge, Sn, and MgO, the method bypasses more than 99% of first-principles calculations, accelerates simulations by a factor of thousand, and quantitatively reproduces first-principles melting points.

  • Problem

    First-principles phase-transition predictions require substantial computational resources, and conventional machine-learning force-field generation needs large datasets and complex trial-and-error selection.

  • Method

    The paper develops an on-the-fly machine-learning force-field method using Bayesian uncertainty estimates to decide when first-principles calculations are needed during molecular dynamics.

  • Results

    Applied to Al, Si, Ge, Sn, and MgO, the method generates accurate force fields and quantitatively reproduces first-principles melting points.

  • Takeaways & Limitations

    The method enables efficient force-field generation for metallic, covalent, and ionic materials while substantially reducing first-principles calculations and simulation cost.

Abstract

from arXiv · show

An efficient and robust on-the-fly machine learning force field method is developed and integrated into an electronic-structure code. This method realizes automatic generation of machine learning force fields on the basis of Bayesian inference during molecular dynamics simulations, where the first principles calculations are only executed, when new configurations out of already sampled datasets appear. The developed method is applied to the calculation of melting points of Al, Si, Ge, Sn and MgO. The applications indicate that more than 99 \% of the first principles calculations are bypassed during the force field generation. This allows the machine to quickly construct first principles datasets over wide phase spaces. Furthermore, with the help of the generated machine learning force fields, simulations are accelerated by a factor of thousand compared with first principles calculations. Accuracies of the melting points calculated by the force fields are examined by thermodynamic perturbation theory, and the examination indicates that the machine learning force fields can quantitatively reproduce the first principles melting points.

I. INTRODUCTION

First-principles melting-point predictions are computationally demanding, while conventional machine-learning force-field generation requires extensive datasets and manual selection. The paper presents an on-the-fly Bayesian scheme that selects new first-principles calculations during molecular dynamics and updates the force field automatically.

  • Motivation: Predicting first-order phase transitions from first-principles calculations is a long-standing condensed-matter problem requiring substantial computational resources.Direct molecular dynamics and Monte Carlo simulations are difficult because phase transitions occur over long timescales.
  • Motivation: Conventional machine-learning force fields often require carefully selected datasets containing 2000–12000 structures and many trial-and-error optimization steps.These requirements have limited applications mainly to a few simple materials.
  • Approach: On-the-fly force-field methods use Bayesian uncertainties to decide whether predicted energies, forces, and stresses are reliable enough to advance molecular dynamics.Large uncertainties trigger first-principles calculations whose results refine the force field; small uncertainties allow the simulation to continue using predictions.
  • Approach: The proposed algorithm generates force fields during molecular dynamics while systematically selecting data across wide phase spaces.The workflow predicts observables, conditionally performs first-principles calculations, updates reference datasets, and advances atomic positions and velocities.
  • Contribution: The method is designed as an efficient and robust algorithm for liquid–solid phase transitions across a wide variety of materials.Its required components were implemented in the Vienna Ab initio Simulation Package.

B. Descriptor

The method represents each atom’s local environment with rotationally invariant descriptors derived from neighboring-atom distributions. Radial and angular information are combined through basis expansions to characterize local structure for force-field construction.

  • Local environment: The total potential energy is decomposed into local atomic energies, each determined by the environment surrounding atom i.The local environment is represented by the probability density of neighboring atoms within a cutoff radius Rcut.
  • Radial and angular descriptors: A radial distribution descriptor captures pairwise distances from the central atom but lacks angular information and can therefore map distinct environments to the same representation.The radial descriptor is illustrated in Fig. 2(a), while the ambiguity arises because different probability densities can yield identical radial distributions.
  • Radial and angular descriptors: Angular information is added through the probability of finding two neighbors at distances r and s separated by angle θ around the central atom.This angular distribution is schematically shown in Fig. 2(b).
  • Basis expansion: The angular distribution function is equivalent to the power spectrum used in practical GAP applications and can be obtained by expanding the neighbor density in radial and angular basis functions.The angular basis uses spherical harmonics, while the radial basis functions satisfy an orthonormal relation.
  • Basis expansion: The resulting coefficients retain the angular-distribution information and provide a descriptor-based representation of the local environment.The text states that the coefficients p^i_nνl contain the same information as the angular distribution defined earlier.

C. Potential energy, Gaussian approximation potential

The potential energy is modeled with a Gaussian approximation potential using kernel similarities between local configurations and reference environments. The implementation combines SOAP-based descriptors with computational optimizations for on-the-fly molecular-dynamics use.

  • Potential energy model: Each local energy is represented as a function of descriptor coefficients, and the total energy is assembled from contributions associated with local configurations.The coefficient vector X^i collects the descriptor coefficients for a specific local configuration.
  • Gaussian approximation potential: The Gaussian approximation potential selects NB local reference structures and represents configurations through corresponding coefficient sets.This approach was adopted because a suitable explicit functional form is generally difficult to determine.
  • Kernel representation: The kernel measures similarity between an environment of interest and reference configurations, approaching unity for similar configurations and decreasing for dissimilar ones.The kernel is constructed from normalized descriptor vectors and weighting parameters.
  • Kernel representation: The kernel contains a pairwise linear term for radial interactions and a nonlinear many-body term identified with SOAP.The radial term is suited to long-range interactions such as Coulomb and Lennard-Jones interactions, while the SOAP term supplies nonlinear many-body interactions.
  • Implementation choices: For all materials in this study, the SOAP kernel is used with β(2)=0, while an application including the radial descriptor is reported elsewhere.This choice was made to enable comparison with previous SOAP-based studies.
  • Implementation choices: Spline interpolation accelerates coefficient and derivative evaluation by tabulating a radial function on a mesh once and reusing its interpolated values during training or molecular dynamics.The implementation also uses normalized spherical Bessel functions as radial basis functions.

D. Fitting of energy, forces and stress tensor and their uncertainty

The method fits energies, forces, and stress tensors with Bayesian linear regression and uses predictive uncertainty plus spilling factors to assess force-field reliability. These estimates guide adaptive sampling during on-the-fly simulations.

  • Fitting: The force field fits dimensionless FP energies per atom, forces, and stress-tensor components as linear functions of descriptor-space coefficients.Reference observations are scaled by their dataset standard deviations before fitting.
  • On-the-fly generation: The on-the-fly scheme performs FP data generation, selection, and parameter optimization during MD rather than preparing broad datasets in advance.Uncertainty evaluation is identified as a key component of the method’s efficiency.
  • Bayesian regression: Bayesian linear regression determines coefficients from limited noisy FP datasets, with regularization parameters optimized through evidence approximation.The posterior maximum coincides with ridge regression, with the noise-to-prior-variance ratio acting as the Tikhonov parameter.
  • Uncertainty estimation: Bayesian inference produces posterior predictions whose mean gives the predicted quantities and whose variances provide prediction uncertainties.The posterior is obtained by propagating the coefficient distribution to the predicted energy, forces, and stress tensor.
  • Spilling-factor assessment: The spilling factor measures the density of nearby local reference configurations and approaches zero for strong overlap or unity when overlap is absent.Unlike Bayesian error, it is more sensitive to structural differences but can suffer numerical instability without regularization.
  • Spilling-factor assessment: Only 2% of FP calculations were triggered by the spilling-factor criterion, which was therefore less relevant than the Bayesian-error criterion in these applications.The spilling-factor threshold was set to 0.02, and it was exceeded only for liquid and interfacial MgO.

E. Decision to perform FP calculation

The decision procedure combines recent sampling history with Bayesian force errors and spilling factors to decide whether a new FP calculation is needed. It skips overly dense sampling and executes FP calculations when estimated errors exceed their criteria.

  • Decision criteria: Within 10 MD steps of the previous sampling event, the procedure skips FP calculations to avoid overly dense sampling in a narrow phase space.After that interval, it evaluates the estimated errors and sampling history.
  • Decision criteria: An FP calculation is performed when the maximum Bayesian force error or spilling factor exceeds its chosen threshold; otherwise, it is skipped.The flowchart uses vectors of per-atom Bayesian force errors and spilling factors together with infinity norms.
  • Threshold adaptation: The spilling-factor threshold is fixed at 0.02, whereas the Bayesian-error threshold is determined automatically during training.The Bayesian threshold is initialized to zero and updated from recent stored maximum errors when their relative variability is sufficiently small.

F. Sparsification and data selection

The data-selection procedure reduces retraining and dataset size by batching FP results, selecting uncertain local configurations, and applying CUR-based sparsification.

  • Retraining: The force field is typically retrained after 5 FP calculations or earlier when estimated errors reach twice their determined criteria.Batching avoids retraining after every individual FP calculation.
  • Local-configuration selection: Local configurations with Bayesian force errors or spilling factors above threshold are selected before a second sparsification step addresses numerical overcompleteness.The procedure then uses a CUR algorithm and leverage scores based on small eigenvalues of the kernel matrix.
  • CUR sparsification: Configurations associated with kernel-matrix eigenvalues below 10^-10 are discarded according to descending leverage scores.The number discarded equals the number of eigenvalues below that threshold.
  • Dataset filtering: Structure datasets contributing no selected local reference configurations are removed to reduce computational cost and memory usage.This final filtering step follows local-configuration sparsification.

A. Training conditions

Training datasets were collected on the fly from solid, liquid, and interfacial systems for Al, Si, Ge, Sn, and MgO. The simulations used material-dependent cells, cutoffs, functionals, and time steps.

  • Systems and datasets: Reference configurations and datasets were collected during MD simulations of solid, liquid, and interfacial systems.Al used 108-atom solid and liquid cells and a 144-atom interface; other materials used 64-atom solid and liquid cells.
  • Electronic-structure settings: Plane-wave cutoff energies were 325, 325, 225, 135, and 520 eV for Al, Si, Ge, Sn, and MgO, respectively.Training used LDA, PBE, PBEsol, and SCAN functionals.
  • Evaluation: Mean absolute errors were evaluated for energy per atom, forces, and stress tensors on 200 solid and liquid configurations at the melting points.The reported units are meV atom^-1, eV Å^-1, and GPa, respectively.

B. Efficiency and accuracy

The on-the-fly force fields substantially reduce first-principles data-generation costs while accelerating molecular dynamics and maintaining accurate energy, force, and stress predictions. Melting-point calculations use interface pinning, with thermodynamic perturbation theory evaluating deviations from first-principles data.

  • Efficiency: More than 99% of first-principles calculations were skipped, reducing force-field-generation computational time by more than a factor of 200.This enabled efficient collection of reference configurations across a wide phase space.
  • Efficiency: Typically fewer than 500 structures and fewer than 1000 local reference configurations were used in the reference datasets.These datasets were smaller than those used in previous studies and were trained for liquids, solids, and interfaces.
  • Efficiency: Force fields accelerated molecular dynamics simulations by factors of 2000 to 5000 compared with first-principles calculations.The acceleration followed efficient sampling of local reference configurations.
  • Accuracy: Mean absolute errors were 5.5 meV atom^-1 for energies, 0.07 eV Å^-1 for forces, and 0.18 GPa for stress tensors.The corresponding root mean square errors were 6.2 meV atom^-1, 0.09 eV Å^-1, and 0.27 GPa, respectively.
  • Melting-point calculations: Interface pinning determines the melting temperature from the point where the chemical-potential difference between solid and liquid becomes zero.A biased constant-temperature, constant-pressure molecular-dynamics simulation constrains an interfacial order parameter between the solid and liquid values.
  • Melting-point calculations: The entropy of fusion and melting-curve slope were evaluated from enthalpy and volume differences using the Clausius–Clapeyron relation.These thermodynamic properties were compared with previously reported and experimental values.
  • Accuracy: Thermodynamic perturbation theory evaluates how deviations between machine-learning and first-principles potentials affect the melting points.First-principles calculations are performed on selected structures from liquid and solid machine-learning-force-field trajectories.

IV. RESULTS: MELTING POINTS

The on-the-fly machine learning force fields reproduce melting points and related thermodynamic properties across Al, Si, Ge, Sn, and MgO, with accuracy depending on the electronic-structure functional. Thermodynamic perturbation corrections generally improve agreement, while SCAN provides the most balanced melting-point accuracy among the tested functionals.

  • Melting-point accuracy: Thermodynamic perturbation corrections improve agreement between MLFF melting temperatures and reported DFT results for several materials and functionals.The corrections are more strongly correlated with total-energy errors than with force or stress-tensor errors.
  • Functional dependence: For Si, LDA underestimates experiment, PBE improves the value but remains low, and SCAN slightly overestimates while agreeing best.PBEsol strongly underestimates the melting point.
  • Functional dependence: For Ge, PBE, PBEsol, and SCAN follow a trend similar to Si, except that SCAN also underestimates the melting point.For Sn, only SCAN yields a reported melting point, and it still underestimates the melting temperature.
  • Thermodynamic properties: The method evaluates melting points, entropies of fusion, volume changes, and melting-curve slopes for the five studied materials.The latter quantities are compared with prior calculations and experimental results.
  • Structural-energy correlation: The melting-point trend across functionals and materials is reasonably correlated with the DFT energy difference between α-tin and β-tin structures.For Sn, the relevance of this energy difference is questioned because melting starts from β-tin, although an approximate linear relation remains.
  • Functional dependence: For Al, LDA, PBEsol, and SCAN closely reproduce experiment, whereas PBE underestimates the melting point.The same functional pattern is reported for MgO.
  • Computational efficiency: The on-the-fly scheme generates accurate MLFFs for metallic, covalent, and ionic materials while accelerating molecular dynamics by more than two orders of magnitude during learning.For large unit cells, the generated force fields are more than three orders of magnitude faster than first-principles calculations.
  • Functional dependence: SCAN consistently outperforms the tested semilocal functionals for melting-temperature predictions and is the most satisfactory overall.The paper judges SCAN to provide the most balanced accuracy across metallic, covalent, and ionic materials.

Appendix A: Radial basis functions

The appendix constructs radial basis functions for expanding the radial components of atomic descriptors. Their parameters and number are selected automatically to reproduce the target radial functions within prescribed accuracy.

  • Radial basis construction: The radial functions f_lm(r, r_ij) arise by expanding broadened atomic distributions in spherical harmonics and radial functions.The angular factors use spherical harmonics evaluated for the neighbor and position directions.
  • Basis-size selection: The number of radial basis functions N_l^R is determined to reproduce f_lm(r, r_ij) within a predefined accuracy.In implementation, radial functions are evaluated on a grid before molecular-dynamics simulations.
  • Radial basis construction: The radial part is expanded in basis functions χ_nl(r) = j_l(q_nl r), with q_nl chosen so j_l(q_nl R_cut) = 0.The basis functions therefore satisfy a boundary condition at the cutoff radius.
  • Basis-size selection: The broadening width σ_atom determines the required number of radial basis functions once the target reproduction accuracy is fixed.The implementation uses an error tolerance of ±0.02 for reproducing the original radial functions.
  • Bayesian formulation: The Bayesian formulation assumes Gaussian deviation and prior distributions for deriving posterior distributions over model coefficients and predictions.The likelihood p(Y|w) represents the probability of observing first-principles data Y given coefficients w.
  • Bayesian formulation: Completing the square converts the coefficient posterior into a Gaussian form, from which the predictive posterior p(y|Y) is obtained.The predictive distribution is likewise represented as a Gaussian distribution.

Appendix C: Maximization of evidence function

The appendix describes evidence-function maximization and numerical stabilization, then explains a modified CUR procedure for selecting local reference configurations. The modification targets configurations associated with small eigenvalues rather than retaining only dominant components.

  • Evidence-function maximization: The evidence function is maximized with respect to σ2_w by simultaneously solving equations derived from ∂p/∂σ2_w = 0.
  • Numerical stabilization: Eigenvalues smaller than 10^-10 are excluded to avoid numerical instabilities caused by a potentially non-positive-definite matrix.The instability arises when all eigenvalues are used in the computations.
  • CUR formulation: The matrix K is formed from kernel elements Kij for all candidate local reference configurations and is then diagonalized.Both matrix dimensions equal the number of candidates.
  • Modified CUR selection: The modified CUR algorithm discards columns strongly correlated with the N_low eigenvectors having small eigenvalues lξ.The corresponding local configurations are disregarded because sparsification usually discards only a small number of configurations.
  • Modified CUR selection: Statistical leverage scoring measures each column’s correlation with the selected eigenvectors.
Loading 1904.12961v3…