Sparse Gaussian process based machine learning first principles potentials for materials simulations: Application to batteries, solar cells, catalysts, and macromolecular systems

Willow, Soohaeng Yoo; Hajibabaei, Amir; Ha, Miran; Yang, David ChangMo; Myung, Chang Woo; Min, Seung Kyu; Lee, Geunsik; Kim, Kwang S.

doi:10.1063/5.0231265

To design new materials and understand their novel phenomena, it is imperative to predict the structure and properties of materials that often rely on first-principles theory. However, such methods are computationally demanding and limited to small systems. This topical review investigates machine learning (ML) approaches, specifically non-parametric sparse Gaussian process regression (SGPR), to model the potential energy surface (PES) of materials, while starting from the basics of ML methods for a comprehensive review. SGPR can efficiently represent PES with minimal ab initio data, significantly reducing the computational costs by bypassing the need for inverting massive covariance matrices. SGPR rank reduction accelerates density functional theory calculations by orders of magnitude, enabling accelerated simulations. An optimal adaptive sampling algorithm is utilized for on-the-fly regression with molecular dynamics, extending to interatomic potentials through scalable SGPR formalism. Through merging quantum mechanics with ML methods, the universal first-principles SGPR-based ML potential can create a digital-twin capable of predicting phenomena arising from static and dynamic changes as well as inherent and collective characteristics of materials. These techniques have been applied successfully to materials such as solid electrolytes, lithium-ion batteries, electrocatalysts, solar cells, and macromolecular systems, reproducing their structures, energetics, dynamics, properties, phase-changes, materials performance, and device efficiency. This review discusses the built-in library universal first-principles SGPR-based ML potential, showcasing its applications and successes, offering insights into the development of future ML potentials and their applications in advanced materials, catering to both educational and expert readers.

I. INTRODUCTION

A breakthrough in functional materials is driving innovation across various technological fields, spanning from clean energy to information processing. Historically, the discovery of new materials has been hindered by costly trial-and-error methods. Meanwhile, advancements in deep learning, particularly in language, vision, and biology, have demonstrated remarkable predictive capabilities with the accumulation of data and computational power. Leveraging these trends, researchers have developed various machine learning (ML) methods. In particular, in developing new functional materials, the prediction of the structure and properties of the desired materials is possible from various ML approaches to possible materials candidates. One noteworthy achievement is the recent success of ML-based design in creating new solid-state materials.^1,2

To achieve accurate predictions of diverse material systems, conducting simulations using ab initio calculations such as first-principles theory is the most desirable. Consequently, the most comprehensive approach to designing new functional materials involves employing simulation methods that range from simple to complex material systems, operating at the quantum mechanical level. However, ab initio calculations demand substantial computational resources and are often constrained to systems comprising only few hundred atoms. To overcome these limitations, numerous efforts have been directed toward utilizing ML approaches to generate first-principles level interatomic potentials (IAPs) and forces for materials. These endeavors span a wide array of ML approaches to generate the first-principles level IAPs of large-scale materials.^3–35

Diverse ML methods allow us to extend the high quality of ab initio level calculations to large systems and long time periods. Indeed, density functional theory (DFT) combined with ML has enabled fast and realistic computer experiments, which correlate very well with the real world. Notable methods are neural network (NN) representations,^36–42 Bayesian linear regression (BLR)^43–46 and Gaussian process regression (GPR),^25,47–49 gradient-domain⁵⁰ and symmetrized gradient-domain ML,⁵¹ deep potential,⁵² compressed sensing,⁵³ sparse Gaussian process regression (SGPR),⁵⁴ etc.

Based on these various ML methods, there have been diverse applications for materials simulations such as inorganic crystals,^55–57 metal particles,^58,59 solid electrolytes,^54,60–63 batteries,^64–66 electrocatalysts,^67–70 perovskite solar cells,^71–73 molecular liquids,^74–78 organic/bioorganic molecular systems,^79–81 macromolecules,⁸² supramolecules, and polymers.⁸³ In particular, the SGPR rank reduction^84,85 was able to accelerate a DFT calculation of ∼1000 atoms by a factor of 10⁴, resulting in accelerated simulations of silicon crystals, sulfide solid electrolytes as well as others.⁵⁴

In this context, this topical review outlines ML representations of the potential energy surface (PES) of materials using non-parametric Bayesian and GPR algorithms, which have the advantage of estimating the possible errors in the predicted potential and force fields. These representations can effectively capture the PES with a minimal amount of ab initio data, contrasting with NNs that typically require extensive datasets for training their hyper-parameters. The computational efficiency of the GPR algorithm is primarily determined by the inversion of the covariance matrix between data components, such as energy and forces. Addressing this computational bottleneck, SGPR representations have been developed for PES, wherein large datasets are projected into smaller inducing sets, bypassing the need for inverting massive covariance matrices and significantly reducing computational costs. Several orders of magnitude improvement in computational cost is achieved. Moreover, an optimal adaptive sampling algorithm has been devised to facilitate on-the-fly data generation for regression, in conjunction with molecular dynamics (MD) simulations. Through hierarchical combinations of expert models, universal potentials have been generated, enabling the modeling of large-scale complexity via a combinatorial approach.⁶³

As a demonstration, this review presents examples such as the reproduction of melting, glass-crystallization, and phase-change for silicon materials and solid electrolytes. The ab initio level MD simulations using on-the-fly ML of IAPs have been conducted to survey Li diffusivity in hundreds of crystals, potentially applicable as electrolytes for all-solid-state batteries. The models generated for these crystals can be amalgamated to create more general and transferable models through the hierarchical combination of expert models, which generate universal potentials, facilitating the simulation of all these materials and their analogs. This leads to applications for various systems, including batteries, solar cells, catalysts, and hydrocarbon macromolecular systems, paving the way for the eventual evolution of ab initio-level force fields for general materials systems. These universal potentials can be useful in designing novel functional molecular/materials systems as well as predicting their structures and properties. The continuous accumulation of these universal potentials could eventually evolve toward super-AI by self-learning process of diverse IAPs, which would lead to the prediction of diverse undiscovered materials. This comprehensive review offers in-depth insights into the evolution of ML potentials and their practical implementations in advanced functional materials. It is crafted to accommodate readers at both introductory and advanced levels of comprehension, ensuring accessibility for a wide range of audiences.

II. THEORETICAL BACKGROUND

A. Introduction to machine learning (ML) methods

Ab initio simulations are an integral part of our understanding of the materials physics and chemistry alongside experiments. The major bottleneck for application of the ab initio simulations is the computational cost, which scales exponentially with the size of system. If $N_{e}$ is the number of electrons, scaling of the computational cost of different ab initio methods varies from $O (N_{e}^{3})$ to higher powers. DFT is one of the most popular and efficient ab initio methods with a nice trade-off between accuracy and computational cost.⁸⁶ The cost for DFT calculations based on the Kohn–Sham scheme, which requires diagonalization to obtain Kohn–Sham orbitals, scales as $O (N_{e}^{3})$ or at best $O (N_{e}^{2} log N_{e})$ unless one uses the fragmental calculations based on divide-and-conquer algorithms.^87,88 From a practical point of view, DFT calculations are affordable only for $\sim 100$ atoms with reasonable computational resources. In many cases for simulating large systems, force fields are used, which express the potential energy in terms of distances, bonds, angles, dihedrals, etc. Force fields are specially popular in simulating biological systems where large macromolecules such as proteins need to be simulated. The parameters in the force field are often set by fitting to limited experimental data or ab initio results. However, these force fields are of varying transferability between different thermodynamic environments and sometimes between the calculation of different properties within the same environment. In many cases, a force field cannot even satisfactorily fit ab initio data because of the many-body quantum effects, which cannot be reduced to bonds, angles, etc.^89,90

As an alternative to empirically fitting the force fields with simple analytical forms, ML techniques have recently emerged as a powerful tool for representing the mathematical entities used in ab initio modeling.¹⁵ Using ML, the potential energy is learned from a set of n ab initio data points. ML reduces the scaling of the computational cost for energy and force calculations to be linear in the number of atoms in the system. This allows us to apply a potential, which may be trained with the data obtained for small systems, to arbitrarily large systems. Currently, the most widely used ML methods are variations of NNs,³⁶ Gaussian approximation potentials (GAPs),⁴⁷ and BLR methods.^43,53,91 In spite of their success stories, each method has shortcomings. NNs and linear methods need large datasets for training. GAPs can be used for robust modeling with small datasets, but their computational cost for training scales as n³ with the size of the data. Here we discuss the SGPR algorithm, which shows very attractive traits in comparison to the previous ML methods. Moreover, we address the issue of optimal data generation, which may be even more important than the regression algorithm. Combined with SGPR, we also introduce an adaptive sampling algorithm that can build robust potentials with a surprisingly small number of DFT data.

The capabilities of this methodology are demonstrated with applications to superionic conducting solid electrolytes. Li diffusivity is studied in dozens of solid electrolytes. For a certain sulfide based electrolyte, the melting, melt-quench, and diffusivity in the glass phase and glass crystallization are simulated, which are found to be consistent with experimental observations. Additionally, since some of the simulated electrolytes had similar chemical elements, the ML models generated for these materials can be combined to give a universal potential for that subset of elements.

In many simulations, even a computationally fast potential is not enough. For instance, slow dynamics or rare events may require inconceivably long molecular dynamics or Monte Carlo simulations in order to achieve statistical certainty. In general, there are two scales that determine the size of the simulation cell and the period for simulation. The size of the simulation cell must be much larger than the length scale of the spatial correlations. Similarly, the period of the simulation must be much longer than the temporal autocorrelation functions. In these situations, enhanced sampling algorithms are necessary in addition to fast potentials. Although the temporal autocorrelation function oscillates long in many cases in chemistry (compared with $\sim 1 ps$ for simple liquids⁹²), the spatial correlation length is often small (compared, for example, with the Debye length $\sim 10 Å$ for salt concentrations in biological fluids⁹³). Two exceptions exist where the spatial correlation length is divergent. The first is at the critical point of continuous phase transitions and the second is when a system exhibits quasi-long-range order. In this regard, fast and accurate ML models for first-principles potential energy are critical for large-scale materials simulations. Thus, in this section, we address the ML representations of the potential energy surface (PES), including NNs, GPR, and BLR, and then discuss SGPR along with the adaptive sampling algorithms.

As for the ML representations of the PES, the inputs are the atomic coordinates and the targets are the potential energy and forces obtained from ab initio calculations such as DFT. A set of (inputs, targets) constitute the data that will be used for modeling. Assembling the data is an important topic by itself. The data largely determine the domain of validity for the generated model and the computational cost. Therefore, the generation of optimal data might be the most important step for wider transferability and acceptance of the simulation results. Here, we will assume that the data are generated beforehand and we wish to model the existing data. We will discuss the general aspects of ML algorithms and briefly introduce the state of the art methods. Optimal generation of the data by algorithms such as on-the-fly sampling or active learning will be discussed later.

B. Notations

Nomenclature will be followed hereafter. It should be noted that symbols like $D (\cdot)$ and $K (\cdot, \cdot)$ are used to distinguish functions and kernels from coefficients, vectors, etc.

C. Representation of the potential energy

Let

R = {{\vec{r}}_{i}}_{i = 1}^{N}

represent the coordinates of N atoms in a physical system. A set of descriptor functions

{D^{α}}

is defined, which maps R into a surrogate space called the descriptor set

{d^{α}}

⁠, where

d^{α} = D^{α} (R)

with a descriptor index α. The essential property of the descriptor functions is that they should be continuous and differentiable with respect to atomic coordinates. Ideally, there should be one-to-one mapping for symmetrically distinct physical systems and their associated descriptor sets. Symmetry in this context refers to operations such as rotations and translations, which leave the potential energy invariant. In ML terminology, the descriptors are often referred to as feature vectors. The potential energy of the system is then expressed as

E = \sum_{α} E^{α} (d^{α}),

(1)

where

{E^{α}}

is the associated set of latent energy functions (LEFs) for the descriptors

D^{α} \equiv E^{α}

⁠. The main goal of ML is to find representations for LEFs using a dataset of ab initio potential energy and forces. For this, LEFs are represented by regression algorithms such as linear regressions, NN, and GPR.

A subcategory of ML representations, which we are most concerned with, is expression of the potential energy in terms of local chemical environments (LCEs), which represents local atomic structures around an atomic position. In this case, it is useful to introduce a cutoff radius r_c as a measure of locality for practical efficiency. The force on atom i depends only on the relative coordinates of other atoms whose distance is less than r_c. This set of neighbors forms the LCE of atom i, formally defined by

ρ_{i} = {{\vec{r}}_{i j}; j \neq i and | {\vec{r}}_{i j} | < r_{c}},

(2)

where

{\vec{r}}_{i j} = {\vec{r}}_{j} - {\vec{r}}_{i}

⁠. Then, the descriptor for atom i depends on its LCE

d_{i} = D (ρ_{i})

and the potential energy is expressed as

E = \sum_{i = 1}^{N} E (d_{i}) .

(3)

Additionally, the force on atom k is

f_{k}^{μ} = - \frac{\partial E}{\partial r_{k}^{μ}} = - \sum_{i = 1}^{N} \frac{\partial E (d_{i})}{\partial r_{k}^{μ}},

(4)

where

r_{k}^{μ}

are the coordinates of atoms k. Derivatives of

E

with respect to atomic coordinates are obtained using the chain rule

\frac{\partial E (d_{i})}{\partial r_{k}^{μ}} = \frac{\partial E (d_{i})}{\partial d_{i}} \frac{\partial d_{i}}{\partial r_{k}^{μ}} .

(5)

In the following, we will assume the representation given by Eq. , but switching to the generalized form given by Eq. is often as simple as reinterpretation of the indices in the given formulas.

The central assumption in the above formulation is that of locality. This assumption is exercised in practice by choosing a finite cutoff radius r_c. There is a trade-off in choosing the cutoff. Choosing a large cutoff dilutes the features of LCEs and increases the computational cost. On the other hand, by choosing a small cutoff, the effect of atoms beyond the cutoff will appear as random forces that cannot be explained. Therefore, the cutoff should be chosen by balancing the errors of the regression algorithm and the random noise due to the atoms beyond the cutoff. Such a principled approach has been taken for ML of the potential energy for amorphous carbon.⁹⁴

D. Complexity of machine learning models

The complexity of a ML model is roughly proportional to the number of parameters present in the model. We are interested in models that can describe the training data with sufficient accuracy, but the main goal is to use the model for reliable predictions of the latent function at unseen inputs. If this is the case for a model, it is said that the model generalizes well. Otherwise, if a model describes the training data well but fails to predict the latent function at unseen inputs, then it is said that the model over-fits the training data. Therefore, available data are often split into training and testing sets. Only the training set is used for optimizing the parameters of the model. The testing set is only used to estimate the generalization error. Over-fitting can be viewed as a model memorizing the training data instead of finding patterns and relationships among the data. As a rule of thumb, a model with more complexity is more likely to over-fit. Therefore, systematic methods exist to reduce the number of parameters used in the model. This task is called regularization and is often carried out by defining loss functions that aim at minimizing or nullifying some of the parameters in the model. If

w = [w_{1}, \dots]

are the parameters/weights in a model, the most common loss functions for regularization are

L_{1} = | | w | |_{1} = \sum_{j} | w_{j} |

(6)

and

L_{2} = | | w | |_{2}^{2} = \sum_{j} | w_{j} |^{2} .

(7)

For the definition of the norms, see Nomenclature The important distinction between L₁ and L₂ regularization losses is that L₁ encourages some of the weights to be exactly zero, even if some other weights remain large, while L₂ shrinks all of the weights to small but finite values. Therefore, L₁ can be included to nullify some of the weights. L₁ is called the least absolute shrinkage and selection operator (LASSO). L₂ is called ridge regression or Tikhonov regularization. Mixing of L₁ and L₂ is the topic for elastic net regression.

E. Neural networks

A NN consists of stacked layers of weights and biases for the propagation of a descriptor vector by a sequence of transformations. Algebraically,

x^{i + 1} = a^{i + 1} (w^{i + 1, i} x^{i} + b^{i + 1}),

(8)

where w and b are the weights and biases, respectively. The first input is the descriptor vector

x^{1} = d

and the final output is a scalar value for the latent function

E

⁠.

a (\cdot)

is called the activation function, which is just a non-linear transformation to enable the NN to model non-linear data. The layers between the input and the output are called hidden layers (Fig. 1). NNs are used extensively for representing the PES.^36,37,52 Here we will briefly discuss the main characteristics of NNs.

FIG. 1.

View large Download slide

Neural network representation of the latent energy function $E$ ⁠.

For training, the parameters (i.e., weights and biases) in a NN are optimized using a flavor of the gradient descent algorithm. A loss function is defined whose minimization yields the optimal parameters. A common loss function is the squared error

L_{E} = \sum_{i = 1}^{n} {(E_{i} - {\tilde{E}}_{i})}^{2},

(9)

where

{E_{i}}_{i = 1}^{n}

and

{{\tilde{E}}_{i}}_{i = 1}^{n}

are the energies predicted by the NN and energies in the training data, respectively. A loss function L_f can be similarly defined with the force data. An important issue to keep in mind is that NNs are massively parameterized; therefore, they can easily over-fit the training data. For regularization, often L₁ or L₂ loss functions are included in the total loss function alongside L_E and L_f.

NNs are specially efficient for big data since the computational cost for training a NN for data of size n scales as $O (n)$ ⁠. On the other hand, it is not easy to train a NN model using a small training set without over-fitting. Therefore, NNs are very useful when large datasets are already available. However, in the context of ML of the PES, it should be noted that the data are obtained from expensive ab initio calculations and are often scarce.

ML potentials based on NNs have greatly improved due to advancements in graph-based algorithms.^{1,17,40,95–100} When the predicted energies are constrained to be invariant under rotation, translation, and permutation, and the predicted forces are constrained to be equivariant under rotation as the NN parameters are updated, the overfitting issues are significantly reduced.

Graph neural networks (GNNs) are useful for modeling systems where data can naturally be represented as a graph, such as atomic systems in materials and molecular simulations. In these models, nodes represent atoms, and edges represent interactions like bonds or forces between them. Two important categories of GNNs in this domain are invariant models (M3GNet,⁹⁶ CHGNet⁹⁷) and equivariant models (NequIP,⁴⁰ MACE¹⁷). Invariant models ensure the output properties are invariant to the rotation, translation, and reflection of the input system, while equivariant models enforce the predictions change in a way that follows the symmetry transformations (like rotations) of the input system. Invariant models are typically simpler and faster, suitable for scalar outputs like total energy, while equivariant models that incorporate geometric symmetry are ideal for tasks requiring directional dependencies (such as forces and torques).

F. Gaussian processes

A Gaussian process is a type of stochastic process that is generalization of the multivariate normal (MVN) distribution to the continuous index set. A stochastic process is a set of random variables

{f (x) : x \in X}

⁠, where

X

is an index set (for instance,

X = R^{d}

⁠, where d represents a dimension). A Gaussian process is a stochastic process such that any finite subset of random variables has a MVN distribution. A MVN distribution with the mean vector m and the covariance matrix K is defined by

p (f) = \frac{1}{\sqrt{{(2 π)}^{n} | K |}} exp (- \frac{1}{2} {(f - m)}^{T} K^{- 1} (f - m)),

(10)

where

f = {[f_{1}, \dots, f_{n}]}^{T}

⁠. Symbolically,

f \sim N (m, K) .

(11)

A Gaussian process with the mean function

M (\cdot)

and the covariance kernel

K (\cdot, \cdot)

is defined such that for any finite set of elements

x_{1}, \dots, x_{n} \in X

⁠, the associated finite set of random variables

f (x_{1}), \dots, f (x_{n})

have a MVN distribution

[\begin{matrix} f (x_{1}) \\ ⋮ \\ f (x_{n}) \end{matrix}] = N ([\begin{matrix} M (x_{1}) \\ ⋮ \\ M (x_{n}) \end{matrix}], [\begin{matrix} K (x_{1}, x_{1}) & \dots & K (x_{1}, x_{n}) \\ ⋮ & ⋱ & ⋮ \\ K (x_{n}, x_{1}) & \dots & K (x_{n}, x_{n}) \end{matrix}])

(12)

or symbolically

f \sim G P (M, K) .

(13)

The covariance kernel

K (x, x')

is a positive definite function, given by

K (x, x') \equiv 〈 f (x) f (x') 〉,

(14)

which results in a symmetric positive definite covariance matrix K such that

K_{i j} = K (x_{i}, x_{j}) .

(15)

Positive definiteness ensures that all of the eigenvalues of the resulting covariance matrix K are positive. For intuition, let

f (x) = \sum_{j} ϵ_{j} ϕ_{j} (x),

(16)

where

Φ (x) = {ϕ_{j} (x)}

is the mapping to decoupled components such that ϵ_j are independent random variables

〈 ϵ_{i} ϵ_{j} 〉 = σ_{j}^{2} δ_{i j}

⁠. Then for the covariance of f at x and

x'

⁠, we have

〈 f (x) f (x') 〉 = \sum_{j} σ_{j}^{2} ϕ_{j} (x) ϕ_{j} (x') .

(17)

The reverse argument is Mercer's theorem, which states that a given kernel can be expressed in terms of its eigenvalues and eigenfunctions (⁠

\int K (x, x') ϕ (x) d μ (x) = λ ϕ (x')

⁠, where

d μ (x)

may be dx or

p (x) d x

⁠)

K (x, x') = \sum_{i} λ_{i} ϕ_{i} (x) ϕ_{i} (x') .

(18)

As an example for the covariance kernel, the radial basis function (RBF) kernel is defined by

K (x, x') = exp (- \frac{| x - x' |^{2}}{2 ξ^{2}}),

(19)

where ξ is a hyper-parameter.

For representing the PES, the latent energy function

E

is taken as the Gaussian process. A detailed formulation of the regression algorithm will be discussed in a later section. Before proceeding, it is useful to also consider the covariance between derivatives of the latent function since we need to consider forces in addition to energies:

〈 \frac{\partial f (x)}{\partial x} f (x') 〉 = \frac{\partial K (x, x')}{\partial x},

(20)

〈 f (x) \frac{\partial f (x')}{\partial x'} 〉 = \frac{\partial K (x, x')}{\partial x'},

(21)

〈 \frac{\partial f (x)}{\partial x} \frac{\partial f (x')}{\partial x'} 〉 = \frac{\partial^{2} K (x, x')}{\partial x \partial x'} .

(22)

1. Multivariate normal distributions

Here we review some of the most important properties of MVN distributions defined by Eq. . Let f and

f'

be random vectors with MVN distributions defined as

f \sim N (m, K),

(23)

f' \sim N (m', K') .

(24)

If f and

f'

have the same size, then

f + f'

is also a MVN distribution given by

f + f' \sim N (m + m', K + K') .

(25)

Let f be partitioned into two parts

f_{A}, f_{B}

such that

f = [\begin{matrix} f_{A} \\ f_{B} \end{matrix}], m = [\begin{matrix} m_{A} \\ m_{B} \end{matrix}], K = [\begin{matrix} K_{A A} & K_{A B} \\ K_{B A} & K_{B B} \end{matrix}], K_{A B} = K_{B A}^{T} .

(26)

Then, the marginal density is defined by

p (f_{A}) = \int_{f_{B}} p (f_{A}, f_{B}) d f_{B} .

(27)

One can show that

p (f_{A})

is a MVN distribution given by

f_{A} \sim N (m_{A}, K_{A A}) .

(28)

Similarly, the conditional density defined by

p (f_{A} | f_{B}) = \frac{p (f_{A}, f_{B})}{p (f_{B})}

(29)

is a MVN distribution (see the Appendix for the derivation) given by

f_{A} | f_{B} \sim N (m_{A} + K_{A B} K_{B B}^{- 1} (f_{B} - m_{B}), K_{A A} - K_{A B} K_{B B}^{- 1} K_{B A}) .

(30)

Since the covariance matrix is symmetric positive definite, it can be decomposed as

K = L L^{T}

⁠, where L is a lower triangular matrix. This factorization is called Cholesky decomposition. Any symmetric matrix has the eigen-decomposition defined as

K = Q Λ Q^{T}

⁠, where columns of Q are eigenvectors of K and Λ is a diagonal matrix whose entries are eigenvalues of K. The inverse of covariance matrix can be expressed as

K^{- 1} = Q Λ^{- 1} Q^{T}

⁠. Then the vector

α = Q^{T} (f - m)

has a Gaussian product distribution

p (α) = \prod_{i} \frac{1}{{(2 π Λ_{i})}^{1 / 2}} exp (- \frac{α_{i}^{2}}{2 Λ_{i}}) .

(31)

2. Gaussian process regression

In GPR, for a given data

{(x_{i}, y_{i})}_{i = 1}^{l}

⁠, we assume

y_{i} = f (x_{i}) + ϵ_{i},

(32)

where

ϵ \sim N (0, σ^{2})

is a random noise, x_i is a columnar (feature) vector of size d, y_i is a scalar (label), and f is the latent function. For briefness,

x = {[\begin{matrix} \underset{____}{} {(x_{1})}^{T} \underset{____}{} \\ ⋮ \\ \underset{____}{} {(x_{l})}^{T} \underset{____}{} \end{matrix}]}_{l \times d}, f = [\begin{matrix} f (x_{1}) \\ ⋮ \\ f (x_{l}) \end{matrix}], e = [\begin{matrix} ϵ_{1} \\ ⋮ \\ ϵ_{l} \end{matrix}], y = [\begin{matrix} y_{1} \\ ⋮ \\ y_{l} \end{matrix}] .

(33)

Considering a valid covariance kernel

K (\cdot, \cdot)

⁠, resulting in the covariance matrix

K_{i j} = K (x_{i}, x_{j})

⁠, we assume a Gaussian process prior for f as

f = G P (0, K) \to f \sim N (0, K) .

(34)

Then, a MVN distribution for the labels will be assumed

y | x = f + e | x \sim N (0, K + σ^{2} I) .

(35)

We can split the data into training and test sets as follows:

y = [\begin{matrix} y_{n} \\ y_{*} \end{matrix}], x = [\begin{matrix} \underset{__}{} x_{n} \underset{__}{} \\ \underset{__}{} x_{*} \underset{__}{} \end{matrix}], K = [\begin{matrix} K_{n n} & K_{n *} \\ K_{* n} & K_{* *} \end{matrix}] .

(36)

Then, Eq. becomes

[\begin{matrix} y_{n} \\ y_{*} \end{matrix}] | x_{n}, x_{*} \sim N (0, [\begin{matrix} K_{n n} + σ^{2} I & K_{n *} \\ K_{* n} & K_{* *} + σ^{2} I \end{matrix}]) .

(37)

By removing the noise from the testing data, we get

[\begin{matrix} y_{n} \\ f_{*} \end{matrix}] | x_{n}, x_{*} \sim N (0, [\begin{matrix} K_{n n} + σ^{2} I & K_{n *} \\ K_{* n} & K_{* *} \end{matrix}]) .

(38)

Then, using the conditional density formula [Eq. ], we can find the predictive distribution for

f_{*}

as

f_{*} | y_{n}, x_{n}, x_{*} \sim N (μ^{*}, ν^{*}),

(39)

where the predictive mean and variance are given by

\begin{matrix} μ^{*} = K_{* n} {(K_{n n} + σ^{2} I)}^{- 1} y_{n}, \\ ν^{*} = K_{* *} - K_{* n} {(K_{n n} + σ^{2} I)}^{- 1} K_{n *} . \end{matrix}

(40)

In contrast to NNs, which are massively parameterized, Gaussian processes are (almost) non-parametric. Often, only a few hyper-parameters are included in the covariance kernel, which need to be optimized. Therefore, GPR works well with small datasets and is almost prone to over-fitting. The schematic for GPR, covariance, and kernel is shown in Fig. 2.

FIG. 2.

View large Download slide

Schematic features of Gaussian process regression and of covariance and kernel.

Meanwhile, for GPR, we need to invert the covariance matrix on the training data, which is computationally costly. If the size of covariance matrix is n, the computational cost for matrix inversion is $O (n^{3})$ ⁠. For reducing this unfavorable scaling of the cost, several approximation methods have emerged that reduce this scaling to $O (n)$ ⁠. We will discuss one of them in a later section.

3. Hyper-parameters

Any selected covariance kernel

K (\cdot, \cdot)

has a set of hyper-parameters ϑ. In training, we tune these parameters to maximize the probability of training labels (marginal likelihood)

y | x, ϑ \sim N (0, K + σ^{2} I) .

(41)

In practice, we maximize logarithm of this probability (log-marginal-likelihood)

log p (y | x, ϑ) = - \frac{1}{2} [y^{T} Γ^{- 1} y + log | Γ | + n log (2 π)],

(42)

where

Γ = K + σ^{2} I

⁠. For gradient ascent algorithms, the partial derivatives are

\frac{\partial}{\partial ϑ_{j}} log p = \frac{1}{2} Tr {(α α^{T} - Γ^{- 1}) \frac{\partial Γ}{\partial ϑ_{j}}},

(43)

where

α = Γ^{- 1} y

⁠.

If the inverse of the matrix is not explicitly required in the calculation of quantities such as $Γ^{- 1} y$ ⁠, we can bypass the matrix inversion and calculate the result α as the solution of linear system $Γ α = y$ ⁠. If Γ is symmetric and positive definite, the Cholesky decomposition can be used to significantly reduce the computational complexity. An important point is worth noting here about $σ^{2}$ ⁠, which is added to the diagonal elements of the covariance matrix. This should not be set to exactly zero, even if our data are noiseless. A small non-zero σ is needed to ensure the numerical stability of the Cholesky factorization and matrix inversion. The reason is that although the covariance matrix is in principle positive definite, an eigenvalue can become 0 due to numerical errors. This makes the determinant 0 and the inverse non-existent.

4. Covariance matrix for energies and forces

For GPR of the potential energy defined by Eq. , the latent energy function

E (\cdot)

is assumed to be a Gaussian process with the covariance kernel

K (\cdot, \cdot)

K (d, d') \equiv 〈 E (d) E (d') 〉 .

(44)

Explicit forms for

K (\cdot, \cdot)

will be discussed later. For two sets of atomic positions R and

R'

⁠, the covariance between respective potential energies is obtained by

α (R, R') = 〈 E (R) E (R') 〉 = \sum_{i, j} 〈 E (d_{i}) E (d'_{j}) 〉 = \sum_{i, j} K (d_{i}, d'_{j}),

(45)

where d_i and

d'_{j}

are the descriptors for the LCEs in R and

R'

⁠. The covariance between a force component and a total energy is expressed as

〈 f_{k}^{μ} (R) E (R') 〉 = - 〈 \frac{\partial E (R)}{\partial r_{k}^{μ}} E (R') 〉 = - \frac{\partial α (R, R')}{\partial r_{k}^{μ}} = - \sum_{i, j} \frac{\partial d_{i}}{\partial r_{k}^{μ}} \frac{\partial K (d_{i}, d'_{j})}{\partial d_{i}},

(46)

where

r_{k}^{μ}

are the coordinates of the atom k in R and the chain rule is applied at the last step. Similarly for two force components, we can write

〈 f_{k}^{μ} (R) f_{l}^{ν} (R') 〉 = 〈 \frac{\partial E (R)}{\partial r_{k}^{μ}} \frac{\partial E (R')}{\partial r_{l}^{ν}} 〉 = \frac{\partial^{2} α (R, R')}{\partial r_{k}^{μ} \partial r_{l}^{ν}} = \sum_{i, j} \frac{\partial d_{i}}{\partial r_{k}^{μ}} \frac{\partial^{2} K (d_{i}, d'_{j})}{\partial d_{i} \partial d'_{j}} \frac{\partial d'_{j}}{\partial r_{l}^{ν}} .

(47)

Let the number of atoms in R and

R'

be N and

N'

⁠, respectively. Then the computational cost for Eq. is

C_{1} \propto N N' .

For Eq. , since

\partial d_{i} / \partial r_{k}^{μ}

is non-zero only if i and k are neighbors, the computational cost is

C_{2} \propto n_{c} N',

where n_c is the number of atoms in the neighborhood of atoms k within the cutoff radius r_c. Similarly, for Eq. , the cost becomes

C_{3} \propto n_{c}^{2} .

The total cost for all covariance calculations, considering that there are 3N and

3 N'

force components in R and

R'

⁠, becomes

C \propto {(1 + 3 n_{c})}^{2} N N' .

(48)

G. Potential energy inference

Here we discuss the algebra for GPR of the ab initio potential energy and forces using a dataset of such calculations. For two sets of atomic positions R and

R'

⁠, we can calculate the covariance matrix as

K_{R R'} = [\begin{matrix} α (R, R') & \dots & - \frac{\partial α (R, R')}{\partial r_{l}^{ν}} & \dots \\ ⋮ & ⋮ \\ - \frac{\partial α (R, R')}{\partial r_{k}^{μ}} & \dots & \frac{\partial^{2} α (R, R')}{\partial r_{k}^{μ} \partial r_{l}^{ν}} & \dots \\ ⋮ & ⋮ \end{matrix}],

(49)

where

α (R, R')

given by Eq. and

r_{k}^{μ}, r_{k}^{ν}

are the coordinates of atoms in R and

R'

⁠. If the numbers of atoms in R and

R'

is N and

N'

⁠, the shape of

K_{R R'}

is

(1 + 3 N) \times (1 + 3 N')

⁠. In addition, let

Y_{R} = [\begin{matrix} {\tilde{E}}_{R} \\ ⋮ \\ {\tilde{f}}_{k}^{μ} \\ ⋮ \end{matrix}]

(50)

be a vector of the length

(1 + 3 N)

containing the potential energy and forces obtained from ab initio calculations for R.

With above definitions, let

R_{1}, \dots, R_{n}

be n distinct ab initio calculation. Using these dataset, we wish to predict the ab initio potential energy and forces for

R_{*}

⁠. Using Eq. , the predictive mean is

Y_{R_{*}} = K_{* n} w,

(51)

w = {(K_{n n} + σ^{2} I)}^{- 1} y_{n},

(52)

where

\begin{matrix} K_{* n} = K_{n *}^{T}, K_{n *} = [\begin{matrix} K_{R_{1} R_{*}} \\ ⋮ \\ K_{R_{n} R_{*}} \end{matrix}], K_{n n} = [\begin{matrix} K_{R_{1} R_{1}} & \dots & K_{R_{n} R_{1}} \\ ⋮ & ⋮ \\ K_{R_{1} R_{n}} & \dots & K_{R_{n} R_{n}} \end{matrix}], y_{n} = [\begin{matrix} Y_{R_{1}} \\ ⋮ \\ Y_{R_{n}} \end{matrix}] . \end{matrix}

(53)

The predictive variance becomes

ν^{*} = K_{* *} - K_{* n} {(K_{n n} + σ^{2} I)}^{- 1} K_{* n}^{T} .

(54)

For improving the inference, we need to maximize

F = - \frac{1}{2} (y_{n}^{T} {(K_{n n} + σ^{2} I)}^{- 1} y_{n} + log | K_{n n} + σ^{2} I |)

(55)

with respect to the hyper-parameters in the covariance kernel

K (\cdot, \cdot)

and possibly the descriptor function

D (\cdot)

⁠.

If all of the systems in the dataset have the same number of atoms N, the size of the covariance matrix that needs to be inverted is

(1 + 3 N) n

⁠. The computational cost for calculating w is

C_{w} \propto {(1 + 3 N)}^{3} n^{3} .

(56)

Using Eq. , the computational cost for construction of the covariance matrix k_nn is

C_{cov} \propto {(1 + 3 n_{c})}^{2} N^{2} n^{2} .

(57)

With

N \sim 100

atoms and reasonable computational resources, in practice only a few hundred ab initio data can be used for full-rank GPR with an acceptable efficiency. Therefore, approximations for reduction of the cost scaling are crucial. Variations of this formalism are used for regression of the ab initio potential energy.^{44,47,49–51,101}

H. Bayesian linear regression

Let

z = {χ_{j}}_{j = 1}^{m}

be a set of descriptors. This set can be a subset of all of the LCEs in n ab initio data

R_{1}, \dots, R_{n}

⁠. For a given descriptor d, the latent function is expressed as

E (d) = \sum_{j = 1}^{m} w_{j} K (d, χ_{j}),

(58)

where

w = {[w_{1}, \dots, w_{m}]}^{T}

are the weights for the descriptors in z and

K (\cdot, \cdot)

is a covariance kernel similar to GPR. For atomic positions given by R, we define

β (R, z) = \sum_{i = 1}^{N} [K (d_{i}, χ_{1}), \dots, = K (d_{i}, χ_{m})],

(59)

where d_i are the descriptors in R, N is the number of atoms, and

β (R, z)

is a

1 \times m

matrix. The potential energy for R becomes

E (R) = \sum_{i = 1}^{N} E (d_{i}) = β (R, z) w .

(60)

Then, the forces acting on the kth atom are obtained from

f_{k}^{μ} = - \frac{\partial β (R, z)}{\partial r_{k}^{μ}} w .

(61)

The instantaneous pressure tensor

P^{μ ν}

is defined as

P^{μ ν} = N k_{B} T^{μ ν} / V + W^{μ ν} / V,

(62)

where

T^{μ ν}

is the instantaneous temperature tensor. To predict the pressure tensor, the internal virial tensor

W^{μ ν}

obtained by

\begin{matrix} W^{μ ν} (R, z) = \sum_{i}^{N} f_{i}^{μ} r_{i}^{ν} = \sum_{i}^{N} [- \frac{\partial β (R, z)}{\partial r_{i}^{μ}} w] r_{i}^{ν} . \end{matrix}

(63)

The weights w can be calculated by linear fitting of E(R),

f_{k}^{μ}

⁠, and

W^{μ ν}

to the energies, forces, and virial pressure from ab initio calculations for the data

R_{1}, \dots, R_{n}

⁠. For this, the following matrices are defined:

K_{R m} = [\begin{matrix} β (R, z) \\ ⋮ \\ - \frac{\partial β (R, z)}{\partial r_{k}^{μ}} \\ ⋮ \\ W^{μ ν} (R, z) \\ ⋮ \end{matrix}], K_{n m} = [\begin{matrix} K_{R_{1} m} \\ ⋮ \\ K_{R_{n} m} \end{matrix}] .

(64)

Then, the linear system to be solved is

K_{n m} w = y_{n},

(65)

where

y_{n}

is the vector with the associated ab initio energies, forces, and virial pressures in the data [see Eqs. and ]. The formal solution of this linear system, with a ridge regression regularization, becomes

w = {(σ^{2} I + K_{n m}^{T} K_{n m})}^{- 1} K_{n m}^{T} y_{n},

(66)

where σ is a small regularization parameter [Eq. is the vector that minimizes

| | K_{n m} w - y_{n} | |_{2}^{2} + σ^{2} | | w | |_{2}^{2}

], but computationally it is better to solve Eq. directly. The predictive variance can be calculated from

ν^{*} = σ^{2} I + K_{* m} {(σ^{2} I + K_{n m}^{T} K_{n m})}^{- 1} K_{* m}^{T} .

(67)

The computational cost for w, with the same assumptions that we made in Sec. II G, becomes

C_{w} \propto (1 + 3 N + 9) n m^{2},

(68)

which is significantly smaller than the inversion of the full covariance matrix since typically

m ≪ n

⁠. The cost of calculation of the covariance matrices scales as

C_{cov} \propto (1 + 3 n_{c} + 9) Nnm .

(69)

Determination of the set z is the most challenging step. Naively, the set z can be obtained by collecting all of the LCEs in the data. However, for efficiency, an algorithm can be designed to collect only significantly distinct environments. This method is successfully used for the studies of melting and stability of hybrid metal–organic perovskites.^43,91

I. Compressed sensing

Compressed sensing formalism⁵³ is similar to the linear regression discussed in Sec. II H. The only difference is that instead of using the ridge regression regularization in solving Eq. , the elastic net regression regularization is applied when w is calculated by minimizing the following loss function

L (w) = | | K_{n m} w - y_{n} | |_{2}^{2} + σ^{2} (β | | w | |_{1} + \frac{1 - β}{2} | | w | |_{2}^{2}),

(70)

which reduced to the ridge regression if β = 0. Non-zero β causes some of w_j to vanish. This can be viewed as a systematic method for compressed sampling of the LCE present in the data.

J. Sparse Gaussian process regression (SGPR)

In Secs. II G and II H, we discussed two types of Bayesian regression algorithms: full-rank GPR and BLR methods. In full-rank GPR, the cost for calculation of the covariance matrix scales as $C_{cov} \propto n^{2}$ for a dataset of size n. The cost for evaluation of the regression weights is dominated by inversion of the covariance matrix, which scales as $C_{w} \propto n^{3}$ ⁠. In linear regression methods scaling of both of these costs is reduced to $O (n)$ ⁠. Instead, linear regression methods do not come close in data efficiency, while GPR can represent the potential energy surface with much less data with the same accuracy. In this section, we explore low-rank approximation of the covariance matrix in GPR to reduce the cost scaling to that of the linear regression methods while maintaining the data-efficiency and expressiveness of GPR. This novel algorithm is called SGPR modeling. We also address on-the-fly data sampling methods that are useful to maximize the efficiency of SGPR.

The key for SGPR is low-rank approximation of the covariance matrix in GPR (Figs. 3 and 4). If

z = {χ_{j}}_{j = 1}^{m}

is a reduced set of LCEs, which is sufficient statistics for the LCEs in n ab initio data

R_{1}, \dots, R_{n}

⁠, we can approximate

K_{n n} \approx K_{n m} K_{m m}^{- 1} K_{n m}^{T},

(71)

where K_nn is the full covariance matrix among the data and K_nm is the cross covariance matrix between the data and z. See Eqs. and for explicit forms of these matrices. K_mm is the covariance matrix among z, which is given by

{[K_{m m}]}_{i j} = K (χ_{i}, χ_{j}) .

(72)

FIG. 3.

Schematic of rank reduction of the co-variance matrix of kernel by SGPR ( Knn≈KnmKmm−1KnmT) involved with the inversion of knn, which improves the scalability with data efficient adaptive sampling toward optimal generation of training data and global exploration.

View large Download slide

Schematic of rank reduction of the co-variance matrix of kernel by SGPR (⁠ $K_{n n} \approx K_{n m} K_{m m}^{- 1} K_{n m}^{T}$ ⁠) involved with the inversion of k_nn, which improves the scalability with data efficient adaptive sampling toward optimal generation of training data and global exploration.

FIG. 4.

View large Download slide

Schematic comparison of GPR (left) and SGPR (right). Squares represent atomic configurations. Circles represent compilation of configurations into LCEs around each atom. In GPR, the potential energy for $x^{*}$ is predicted by calculating its covariance matrix with the data X. In SGPR, only the covariance matrix with the inducing LCEs z is needed. Figure taken from Ref. 54. Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

The set z is called an inducing set of LCEs. Similarly, for covariance of the test system

R_{*}

with the data, we approximate

K_{* n} \approx K_{* m} K_{m m}^{- 1} K_{n m}^{T},

(73)

where

K_{* m}

is the covariance matrix between

R_{*}

and z. Using the approximations given by Eqs. and for the predictive mean in Eq. , we can write

Y_{R_{*}} = K_{* m} w,

(74)

w = {(σ^{2} K_{m m} + K_{n m}^{T} K_{n m})}^{- 1} K_{n m}^{T} y_{n},

(75)

where

y_{n}

is the vector of ab initio energies, forces, and virial pressures in the data [see Eqs. and ; refer to Eq. in Appendix]. An interesting observation is that if K_mm in the above equation is replaced by

I

⁠, the predictive mean reduces to that of the BLR method in Eq. . Therefore, except for calculation of K_mm and its inverse, the computational cost is similar to the Bayesian linear regression method. Another point is that if we define

w = K_{m m}^{- 1} e

and substitute into Eq. , it looks like a regression based on

(z, e)

instead of full data. For this reason, the above approximation is also called projected process. The predictive variance becomes

ν^{*} = K_{* *} - K_{* m} K_{m m}^{- 1} K_{* m}^{T} + σ^{2} K_{* m} {(σ^{2} K_{m m} + K_{n m}^{T} K_{n m})}^{- 1} K_{* m}^{T} .

(76)

Finally, w can be viewed as the solution of the following linear system:

[\begin{matrix} K_{n m} \\ σ L^{T} \end{matrix}] w = [\begin{matrix} y_{n} \\ 0 \end{matrix}],

(77)

where L is the Cholesky factor of K_mm (⁠

K_{m m} = L L^{T}

⁠, where L is a lower triangular) and

0

is a columnar vector of zeros with length m.

From another point of view, the SGPR algorithm is just a normal GPR only with the kernel redefined as

K (ρ_{i}, ρ_{j}) \to K_{ρ_{i} m} K_{m m}^{- 1} K_{m ρ_{j}},

(78)

where

K_{ρ_{i} m}

is the covariance matrix between LCE ρ_i and inducing LCE z. That is, the covariance kernel between two LCEs is indirectly calculated using the inducing set z. The covariance loss due to sparsification can be quantified by the spilling factor¹⁰²

s (ρ) = K (ρ, ρ) - K_{ρ m} K_{m m}^{- 1} K_{ρ m}^{T},

(79)

which can also be viewed as the predictive variance of a GPR based on

(z, e)

⁠. Therefore,

s (ρ)

can be utilized in sampling LCEs.

Given that the kernel in SGPR only needs to be calculated between the data and inducing LCEs, a substantial improvement is achieved as compared with full GPR. In a naive comparison, the improvement is $O (n^{2}) \to O (n)$ ⁠, while a more rigorous analysis shows a much better improvement. For example, with 100 configurations (100 atoms each), the improvement is a factor of $O (10^{3})$ ⁠.

K. Variational lower bound for log-likelihood

In full-rank GPR, for optimizing the hyper-parameters, we maximize

log p (y_{n})

⁠, where p is a MVN distribution with the mean 0 and the covariance matrix K_nn. In SGPR, K_nn is approximated by Eq. . In the following, we will refer to the approximated K_nn by

{\tilde{K}}_{n n}

⁠. A straightforward approximation for

p (y_{n})

is to just use

{\tilde{K}}_{n n}

⁠, instead of K_nn. This approximate distribution will be denoted by

\tilde{p}

⁠. Instead of this approximation, Titsias¹⁰³ found a variational lower bound for the true marginal log-likelihood, which is given by

F = log \tilde{p} (y_{n}) - \frac{1}{2 σ^{2}} Tr (K_{n n} - K_{n m} K_{m m}^{- 1} K_{n m}^{T}) .

(80)

Therefore, for optimizing the hyper-parameters in the model, F should be maximized. Adding a LCE to the inducing set z can only increase F. This is very useful for selecting the inducing LCEs. Among several candidates, the one that increases F the most can be selected.

L. Adaptive sampling algorithm

The most important step, regardless of the regression algorithm, is the generation of the data. Early methods^36,47 prepared the data beforehand, which is not efficient. For instance, they used existing ab initio calculations or generated samples by random displacement of atoms. On-the-fly learning methods,^43,44,91,101 estimate the accuracy of the potential and use a criterion for data generation during molecular dynamics simulation, etc. In GPR based algorithms, the predictive variance is a suitable criterion while in NNs the uncertainty can be emulated by training identical NNs with different initial random seeds. Generally, GPR based algorithms are more suited for on-the-fly learning since they have only a few hyper-parameters. However, the size of these potentials and their computational cost grow with time; therefore, the data generated earlier need to be dumped to keep the potential affordable. Thus, for globally accurate potentials, it is critical to devise an optimal adaptive sampling algorithm. We tackle this issue by replacing the predictive variance by another criterion, which we call the geometric criterion. This is motivated from the observation that the predictive variance does not take advantage of the energy and forces data. Assume that a SGPR model is generated using n ab initio data $X = {R_{j}}_{j = 1}^{n}$ and m inducing LCEs $z = {χ_{j}}_{j = 1}^{m}$ ⁠. Aside from hyper-parameter optimization, the potential can be modified in two ways: adding a LCE $χ_{m + 1}$ to the inducing set z and adding a new ab initio data $R_{n + 1}$ to X. Instead of calculating the predictive variance, we calculate the change in the predictive mean of the potential by directly constructing the potential both with and without $R_{n + 1}$ (same for $χ_{n + 1}$ ⁠). If the change in the predictive mean is more than a threshold ϵ, the addition is accepted. Inclusion of the forces data in the regression is critical for this algorithm; otherwise, geometric criterion may become only marginally better than variance criterion.

For adaptive sampling of the data and inducing LCEs, the model is built on-the-fly with MD where the system evolves with SGPR model, and ab initio calculations are carried out actively to correct the model based on the following criterion. The data and unique LCEs in the first step are automatically included in X and z. In the following steps, we try to insert on-the-fly LCEs to z. A LCE $χ_{m + 1}$ is added to z only if $Δ E (χ_{m + 1}) > ϵ$ ⁠, where $Δ E$ is the change in the predictive energy resulting from the inclusion of $χ_{m + 1}$ in z and ϵ is a predetermined threshold. Since it is too costly to try out all LCEs, in the spirit of importance sampling, we try insertion of LCEs based on their spilling factor [see Eq. (79)] in a descending order and terminate at the first unsuccessful insertion trial. If at least one LCE is added to z, it means that the system has potentially crossed into an unfamiliar region. At this stage, in one algorithm (FAST), we calculated the exact ab initio energy and forces and included the snapshot $R_{n + 1}$ in X only if $Δ E (R_{n + 1}) > ϵ$ ⁠. Thus, sometimes the exact ab initio data were regarded as redundant and were rejected by the model. This was devised to keep the number of sampled data (⁠ $\propto$ the computation cost in future steps) as small as possible. In another algorithm (ULTRAFAST), in order to preempt unwanted ab initio calculations, predictions of the model were used as fake ab initio data, and if the insertion was accepted, they were corrected by exact calculations. This algorithm is similar to leave-one-out cross-validation method for active learning.¹⁰⁴

M. AutoForce: A SGPR-based ML package for atomistic simulations

We have developed the Python package AutoForce for generating SGPR models and on-the-fly learning of the PES.¹⁰⁵ This package is coupled with the atomic simulation environment (ASE)¹⁰⁶ and from there it can be coupled to various ab initio software including Vienna Ab initio Simulation Package (VASP),¹⁰⁷ which implements the projector augmented wave (PAW)¹⁰⁸ approach to DFT with Perdew–Burke–Ernzerhof (PBE) and generalized-gradient approximation (GGA) functionals.¹⁰⁹

For reporting the accuracy of ML potentials, we use root mean squared error (RMSE) and coefficient of determination, which is denoted as R² in statistics and is defined by

R^{2} = 1 - \frac{\sum_{i} {(f_{i} - {\tilde{f}}_{i})}^{2}}{\sum_{i} {(f_{i} - \bar{f})}^{2}},

(81)

where

{f_{i}}

are the ab initio forces,

\bar{f}

is their average, and

{{\tilde{f}}_{i}}

are the ML forces. R² indicates the fraction of explained variance.

N. Benchmark test of SGPR-based ML: Silicon

For a benchmark simulation, we created a SGPR potential for bulk Si in that it has been studied by several other ML methods.^{36,44,47,91,101} A cubic cell with 64 atoms was chosen with a $2 \times 2 \times 2$ k-point grid for DFT calculations. Calculations for the remaining part are non-spin-polarized and a kinetic energy cutoff of 500 eV is applied. The potential is generated by on-the-fly sampling with MD simulations at 300 K (for 20 ps) and 1000 K (for 30 ps). Using the algorithm ULTRAFAST, only 13 DFT calculations were performed for data generation and 36 LCEs were sampled as the inducing set. This entire simulation took less than 1 h. A total of 100 snapshots were chosen from trajectories for testing, which resulted in the RMSE = 0.07 eV/AA and $R^{2} = 0.994$ for forces. The phonon spectra with a $5 \times 5 \times 5$ supercell are shown in Fig. 5. This performance is better than all previous ML methods, which used hundreds to thousands of DFT data for generating a potential. The closest competition is on-the-fly active learning,⁴⁴ which reached the RMSE of 0.08 eV/Å with 133 DFT calculations for 10 ps MD at 620 K. The SGPR algorithm, which used much more expensive smooth overlap of atomic positions (SOAP) descriptor and smaller computational resources, was found to be more than ∼300 times faster than the GPR algorithm.

FIG. 5.

View large Download slide

Phonon spectra of bulk Si with a SGPR potential generated using only 13 DFT data and 36 inducing LCEs. Figure taken from Ref. 54. Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

The RMSE of ML potential usually increases with temperature because of larger forces and more disorder. In Si simulation above from 300 to 700 K, the RMSE increases from 0.03 to 0.07 eV/Å. However, at higher temperatures, the larger RMSE can be tolerated because of random forces due to thermal fluctuations. In Ref. 36, a silicon melt simulation with NN errors well above 0.2 eV/Å accurately reproduced the DFT radial distribution function. A minimal requirement is that the distribution of forces error should be Gaussian centered at 0. As such, they act as random forces due to thermal coupling with a heat bath and do not significantly alter the thermodynamics.

O. Universal potentials

When we consider materials with various compositions but with the same elements, each “expert” model can be generated for each composition based on the SGPR algorithm. Then, a general expert model can be obtained for a material system having diverse compositions. As an example, let us consider the Li–S–P materials system (Li_x–S_y–P_z) having four different compositions (⁠ ${Li}_{3} {PS}_{4}, {Li}_{3} {PS}_{6}, {Li}_{48} P_{16} S_{61}, {Li}_{7} P_{3} S_{11}$ ⁠), which are composed of the same set of elements of Li, P, and S. Then, each expert model is generated for each material. By combining these four expert models, one universal expert model can be generated. Three approaches could be utilized:

Product of experts (POE): Average of all of the expert models is calculated for making predictions.
Merging of experts (MOE): The data and inducing LCEs of all expert models are combined by importance sampling using the geometric criterion.
Fusion of experts (FOE): The data and inducing LCEs of all expert models are fully utilized, without importance sampling.

A SGPR model is generated for each material of a certain composition. We may refer to such a model as a local or expert model. An expert model i is obtained from its m_i inducing descriptors and n_i data that are sampled on-the-fly with MD. Then the covariance matrices

k_{m_{i} m_{i}}

and

k_{n_{i} m_{i}}

are built, which yield the weights

w_{i}

needed for predictions [see Eq. ]. Each expert model can be viewed as a block from the universal SGPR model, which is defined by the global covariance matrices

k_{M M}

and

k_{N M}

(with

M = \sum_{i} m_{i}

and

N = \sum_{i} n_{i}

⁠):

k_{M M} = [\begin{matrix} k_{m_{1} m_{1}} & k_{m_{1} m_{2}} & \dots \\ k_{m_{2} m_{1}} & k_{m_{2} m_{2}} & \dots \\ ⋮ & ⋮ & ⋱ \end{matrix}]

(82)

and

k_{N M} = [\begin{matrix} k_{n_{1} m_{1}} & k_{n_{1} m_{2}} & \dots \\ k_{n_{2} m_{1}} & k_{n_{2} m_{2}} & \dots \\ ⋮ & ⋮ & ⋱ \end{matrix}],

(83)

where

k_{m_{i} m_{j}}

and

k_{n_{i} m_{j}}

(for

i \neq j

⁠) are inter-expert cross covariance matrices. Inserting these matrices into Eq. with

Y = {[Y_{1}, Y_{2}, \dots]}^{T}

results in the global weights

W

⁠. Such a model would be able to describe all the materials with the same elements but with different compositions simultaneously, and therefore more transferable. But, building such a universal model is computationally difficult and, even if possible, would become too slow for predictions. The assumption that the off diagonal blocks are zero results in

W = {[w_{1}, w_{2}, \dots]}^{T}

and the prediction of the universal model reduces to simple averaging over expert models. This approximation is referred to as POE. POE becomes reasonable when the expert models are built for completely dissimilar materials. A better approach is to use the expert models as a library. Given a new material x that is not available in the library, we can choose only a subset of expert models, which are similar to x. Algebraically, if

k_{M M}

is built from a subset of experts and

k_{x x} - k_{x M} k_{M M}^{- 1} k_{x M}^{T} \to 0,

(84)

then the chosen subset is a sufficient model. We refer to such a combination as FOE. If such a subset cannot be constructed, then the expert model for x should be added to the library. Since it is possible that similar inducing descriptors are repeated in several expert models, a FOE can be further compressed for increasing the efficiency. This model is referred to as MOE. MOE can be obtained by importance sampling of the FOE data and inducing descriptors using a threshold criterion. Here, the adaptive sampling algorithm described in Ref. 54 is used. In simple terms, POE is obtained by neglecting the off diagonal blocks (cross terms) of the global covariance matrix, MOE by importance sampling of the rows/columns, and FOE by forceful inclusion of all elements.

As an example, here we provide the coefficients of determination R² for the above combinatorial models, i.e., the aforementioned four different crystals of the ternary (Li–P–S) solid electrolyte with the same set of three elements, as listed in Table I. As expected, POE is poorer than expert models. It is included here because in most of the non-Bayesian ML methods, averaging over the expert models is the only available mechanism for combining ML models, which is equivalent to POE in Bayesian methods. The current study provides a quantitative comparison. Thus, although POE can provide a rough ML-IAP without any further ML procedures, it does not generalize well. ML models trained for different domains may not be simply averaged. MOE is as good as the expert models, even though its size (data/inducing) = (102/573) is much smaller than the combined size of the experts (181/1418).⁶³ FOE has the same size as the combined size of the experts but almost in every case is better than the experts. Thus, FOE could be used when affordable; otherwise, MOE offers a scalable solution that can handle more complex systems much faster.

TABLE I.

R² of expert models and their combinations for Li–P–S and Li–Sb–S systems. Adapted from Ref. 63. Reproduced with permission from Hajibabaei and Kim, J. Phys. Chem. Lett. 12, 8115 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

	Li–P–S	Expert	POE	MOE	FOE
A	Li₃PS₄	0.972	0.962	0.973	0.980
B	Li₇PS₆	0.972	0.813	0.966	0.978
C	Li₄₈P₁₆S₆₁	0.953	0.697	0.952	0.969
D	Li₇P₃S₁₁	0.958	0.864	0.955	0.969

	Li–P–S	Expert	POE	MOE	FOE
A	Li₃PS₄	0.972	0.962	0.973	0.980
B	Li₇PS₆	0.972	0.813	0.966	0.978
C	Li₄₈P₁₆S₆₁	0.953	0.697	0.952	0.969
D	Li₇P₃S₁₁	0.958	0.864	0.955	0.969

Similar combinations can also be made for another ternary system such as Li–Ge–S solid electrolyte by accomplishing the same by training a model sequentially in all domains. Then, as a more interesting and useful example, by combining the two expert models for the two ternary systems of Li–P–S and Li–Ge–S, it is possible to generate the expert model for a quaternary system of Li–P–Ge–S solid electrolyte, once the cross term between P and Ge, i.e., the expert model for the P–Ge system is properly provided. If there is no direct contact between these two elements, the cross term could be reliable enough by adding only the long-range interactions between the two elements. As an approximation, such a small effect by the cross term could be even ignored with a certain error range. However, if the two elements are directly in contact, the cross term needs to be taken into account. In the case of MOE, the cross term between the two elements can be more easily and simply calculated than the three elements cases. Overall, the MOE approach can be utilized with a minimal loss of accuracy in energy. Thus, the combinatorial approach utilizing the MOE approach opens the door for strategies like “high level parallelism” or “a library of SGPR potentials,” which can be combined by simple cross covariance calculations to describe increasingly complex materials. Such a combinatorial approach is not trivial by other ML methods.

Indeed, the SGPR models for Li₃PS₄ and Li₄GeS₄ were simply combined and tested on Li₁₀Ge(PS₆)₂ (both configurations in Table II). Due to the presence of different species P and Ge, the combinatorial model cannot be significantly compressed and thus MOE and FOE are equivalent. The combined model of the two ternary crystals described the forces of the quaternary crystal with $R^{2} = 0.961$ ⁠, which is very close to R² of the expert model. This is understandable in view of the similarity of LCEs for the cations P and Ge in these crystals (see Fig. 6). In many promising solid electrolytes, they are present as ${({XS}_{4})}^{- q}$ poly-anionic moieties with X = P, Ge, etc., in the sea of Li⁺ ions and do not directly interact with each other. By combination of the ML models for the simpler crystals, we can simulate the composite ones with a minimal training or even without any additional training.

TABLE II.

22 fastest Li conducting ternary solid electrolytes obtained after filtering of about 300 candidates. Diffusivities are obtained with MLMD in NPT ensemble. The first column is the ID of crystal in the materialsproject¹¹¹ repository. The last column σ_RT is the extrapolated ionic conductivity at 300 K. The quaternary Li₁₀Ge(PS₆)₂ is appended as a special case.⁶³ Adapted from Ref. 63. Reprinted with permission from Hajibabaei and Kim, J. Phys. Chem. Lett. 12, 8115 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

mp-ID	Formula	ML	Diffusion coefficient (Å²/ps) at given temperature (K)								E_a	σ_RT
		R²	400	500	600	700	800	900	1000	1200	(eV)	(mS/cm)
mp-1211324	Li₇PS₆	0.972			0.15		0.29		0.54	0.83	0.18	65.2
mp-1097036	Li₃PS₄	0.972		0.068	0.14	0.28	0.46				0.23	20.2
mp-641703	Li₇P₃S₁₁	0.958			0.12	0.2	0.32		0.61		0.21	20.1
mp-768197	Li₅SbS₄	0.934			0.1	0.18	0.31				0.23	18.4
mp-1147621	Li₃InBr₆	0.924	0.018	0.057	0.13						0.2	16.1
mp-28327	LiGaBr₃	0.903		0.049	0.11						0.21	7.6
mp-756427	Li₅SbS₄	0.933				0.086	0.19	0.26	0.37		0.26	6.0
mp-685863	Li₄SiO₄	0.958					0.093	0.2	0.32	0.53	0.3	3.1
mp-1195718	Li₄SnS₄	0.973					0.18	0.32	0.37	0.74	0.29	2.3
mp-766409	Li₃BiS₃	0.946			0.034	0.077	0.12		0.29		0.27	2.1
mp-29009	Li₂MgBr₄	0.944			0.045	0.096	0.18				0.28	1.4
mp-38684	Li₂MgCl₄	0.958			0.047	0.095	0.21		0.46		0.3	1.3
mp-760415	Li₃SbS₄	0.966					0.2		0.57	0.98	0.33	0.83
mp-1194700	Li₄SnSe₄	0.962					0.15	0.3	0.41	0.8	0.33	0.6
mp-1222582	Li₄GeS₄	0.97					0.14	0.31	0.4	0.93	0.38	0.24
mp-1001069	Li₄₈P₁₆S₆₁	0.953			0.027	0.083	0.17				0.38	0.22
mp-1177520	Li₃SbS₃	0.929			0.018	0.058	0.18		0.39		0.38	0.2
mp-775806	Li₃SbS₃	0.944					0.084	0.16	0.27	0.49	0.36	0.19
mp-38487	Li₇BiO₆	0.987					0.041	0.079	0.12	0.25	0.37	0.17
mp-753720	Li₃BiS₃	0.939					0.08	0.14	0.21	0.51	0.39	0.079
mp-29008	Li₆MgBr₈	0.965			0.019	0.052	0.13	0.23			0.41	0.077
mp-757150	Li₅BiS₄	0.929			0.02	0.12	0.2	0.41			0.46	0.049
mp-696128	Li₁₀Ge(PS₆)₂	0.967			0.16		0.32		0.55	0.76	0.17	84.2
mp-696138	Li₁₀Ge(PS₆)₂	0.97			0.14		0.33		0.55	0.85	0.19	46.6

mp-ID	Formula	ML	Diffusion coefficient (Å²/ps) at given temperature (K)								E_a	σ_RT
		R²	400	500	600	700	800	900	1000	1200	(eV)	(mS/cm)
mp-1211324	Li₇PS₆	0.972			0.15		0.29		0.54	0.83	0.18	65.2
mp-1097036	Li₃PS₄	0.972		0.068	0.14	0.28	0.46				0.23	20.2
mp-641703	Li₇P₃S₁₁	0.958			0.12	0.2	0.32		0.61		0.21	20.1
mp-768197	Li₅SbS₄	0.934			0.1	0.18	0.31				0.23	18.4
mp-1147621	Li₃InBr₆	0.924	0.018	0.057	0.13						0.2	16.1
mp-28327	LiGaBr₃	0.903		0.049	0.11						0.21	7.6
mp-756427	Li₅SbS₄	0.933				0.086	0.19	0.26	0.37		0.26	6.0
mp-685863	Li₄SiO₄	0.958					0.093	0.2	0.32	0.53	0.3	3.1
mp-1195718	Li₄SnS₄	0.973					0.18	0.32	0.37	0.74	0.29	2.3
mp-766409	Li₃BiS₃	0.946			0.034	0.077	0.12		0.29		0.27	2.1
mp-29009	Li₂MgBr₄	0.944			0.045	0.096	0.18				0.28	1.4
mp-38684	Li₂MgCl₄	0.958			0.047	0.095	0.21		0.46		0.3	1.3
mp-760415	Li₃SbS₄	0.966					0.2		0.57	0.98	0.33	0.83
mp-1194700	Li₄SnSe₄	0.962					0.15	0.3	0.41	0.8	0.33	0.6
mp-1222582	Li₄GeS₄	0.97					0.14	0.31	0.4	0.93	0.38	0.24
mp-1001069	Li₄₈P₁₆S₆₁	0.953			0.027	0.083	0.17				0.38	0.22
mp-1177520	Li₃SbS₃	0.929			0.018	0.058	0.18		0.39		0.38	0.2
mp-775806	Li₃SbS₃	0.944					0.084	0.16	0.27	0.49	0.36	0.19
mp-38487	Li₇BiO₆	0.987					0.041	0.079	0.12	0.25	0.37	0.17
mp-753720	Li₃BiS₃	0.939					0.08	0.14	0.21	0.51	0.39	0.079
mp-29008	Li₆MgBr₈	0.965			0.019	0.052	0.13	0.23			0.41	0.077
mp-757150	Li₅BiS₄	0.929			0.02	0.12	0.2	0.41			0.46	0.049
mp-696128	Li₁₀Ge(PS₆)₂	0.967			0.16		0.32		0.55	0.76	0.17	84.2
mp-696138	Li₁₀Ge(PS₆)₂	0.97			0.14		0.33		0.55	0.85	0.19	46.6

FIG. 6.

In all of the above SEs, XS4 poly-anionic moieties are present with X = P or Ge. Figure taken from Ref. 63. Reproduced with permission from Hajibabaei and Kim, J. Phys. Chem. Lett. 12, 8115 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

View large Download slide

In all of the above SEs, XS₄ poly-anionic moieties are present with X = P or Ge. Figure taken from Ref. 63. Reproduced with permission from Hajibabaei and Kim, J. Phys. Chem. Lett. 12, 8115 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

In summary, extensive MD simulations can be carried out with on-the-fly ML for the potential energy using the SGPR algorithm. The models generated with kernel-based methods such as SGPR can be easily combined for more generality and transferability. Therefore, using the previously trained models in the database, it is possible to model new crystals without further training. This approach was demonstrated by creating universal potentials for two ternary systems. Additionally, combinations of two ternary crystals resulted in a realistically accurate potential for a quaternary crystal. With the introduced algorithm, the models generated in this study can be reused, with limited or no training at all, to describe more complex system (such as quaternary electrolytes) by a combinatorial approach.

Recently, universal potentials have been widely developed using graph-based models trained on large databases like the Materials Project.^{1,17,40,95–100} These models aim to cover all elements in the periodic table and handle various configurations, expanding the term “universal potential” to encompass not only elements but also configurational spaces. While the current review discusses mainly elements, nothing in the model formulation limits its applicability to different configurations. The uncertainty threshold of the model can be adjusted to accommodate a broader training configuration space, thereby enhancing its predictive capabilities across diverse materials. The term “universal potential” was first used to refer to extended subsystems within typical machine learning potentials (MLPs), which were tailored for specific materials with fixed compositions.^34,63 Kernel methods are able to address a wide range of material systems, regardless of differences in phase or stoichiometry. These subsystems can be combined into supersystems, enabling them to handle even more diverse materials. A built-in library system, stored in memory, allows for the selection of appropriate kernels for specific systems. Once fully constructed, this library can support various material systems by incorporating cross terms between expert systems, calculated from local structure descriptors. These cross terms are then added to the library, making it more general and versatile over time. Although the compositionally universal potential is currently somewhat limited, the library is continuously expanding and will eventually provide greater generality and applicability for a wide range of materials.

P. Sparse Bayesian committee machine potential

The SGPR algorithm can have difficulties in dealing with large number of inducing sets m as $O (n m^{2})$ in obtaining the solution w of Eq. (74), which involves the inversion of a kernel matrix, with large size of training datasets X and inducing set z of LCEs. In order to mitigate this challenge, Bayesian Committee Machine (BCM)¹¹⁰ is employed to combine SGPR models, each trained independently on different subsets as a local expert (Fig. 7). When building MLP of ice polymorphs, rather than utilizing all ice phases as a training dataset, each ice phase is divided into separate training subsets. Each local expert model is then trained on its corresponding subset. These individually trained local expert SGPR models are combined into a universal BCM MLPs. This approach not only simplifies the training process but also overcomes the computational challenges associated with inverting large kernel matrices, enabling an efficient and accurate universal MLP. The BCM, for instance, combines the predictions from P distinct local expert models, each trained on different data subsets, which reduces the computational scaling to $O (n m^{2} / P^{3})$ for training.

FIG. 7.

View large Download slide

Schematic of active sparse Bayesian committee machine (BCM) potential. (a) The entire dataset is divided into P distinct subsets. Within each subset, a SGPR MLP is trained. Subsequently, a combined version of MLP for the entire dataset is built through the BCM algorithm. The weights $\frac{β_{α}}{σ_{α}^{2}}$ ensure that the predictions from confident models are more outspoken than those from less confident models. (b) Illustration of active BCM algorithm. Each local expert SGPR model undergoes training up to a predefined kernel size. Before surpassing the kernel size limit, our algorithm initiates training for the next SGPR model. This process is sequentially executed, resulting in the assembly of P local SGPR MLPs, BCM(1:P), forming a Bayesian committee. Throughout the MD simulation, this ensemble of P committee MLPs predicts the energy, forces, and virial tensor of the system. Figure taken from Ref. 112. Reproduced with permission from Willow et al., Phys. Chem. Chem. Phys. 26, 22073–22082 (2024). Copyright 2024 Royal Society of Chemistry.

For given subsets of the training dataset

{n_{1}, \dots, n_{P}}

and inducing LCEs

{m_{1}, \dots, m_{P}}

⁠, the BCM¹¹⁰ potential energy is

E = \hat{S} \sum_{α = 1}^{P} \frac{β_{α}}{σ_{α}^{2}} E_{α},

(85)

where the potential energy

E_{α}

of α-th SGPR local expert (⁠

R_{n \in α}

and

z_{α} = {χ_{j}}_{j \in α})

is given by

E_{α} (R) = \sum_{i, (j \in m_{α})} K (ρ_{i}, χ_{j}) w_{j}^{α},

(86)

where

w^{α}

are the training weights of α-th local SGPR model.

The forces

f_{i}^{μ}

and internal virial

W^{μ ν}

of the BCM potential are obtained from

f_{i}^{μ} = \hat{S} \sum_{α = 1}^{P} \frac{β_{α}}{σ_{α}^{2}} f_{i, α}^{μ},

(87)

W^{μ ν} = \hat{S} \sum_{α = 1}^{P} \frac{β_{α}}{σ_{α}^{2}} W_{α}^{μ ν},

(88)

where

f_{i, α}^{μ}

and

W_{α}^{μ ν}

represent the forces and the internal virial of the α-th SGPR local expert, respectively, which are given by

\begin{matrix} f_{i, α}^{μ} = \sum_{j \in m_{α}} {\dot{K}}_{i}^{μ} (R, χ_{j}) w_{j}^{α} and W_{α}^{μ ν} = \sum_{i}^{N} f_{i, α}^{μ} r_{i}^{ν} . \end{matrix}

The maximum of the covariance loss for the α-th subset LCEs, denoted as

σ_{α}^{2}

⁠, is used to weight the α-th committee prediction. The

σ_{α}^{2}

is estimated by

σ_{α}^{2} = max (1 - K_{* m} K_{m m}^{- 1} K_{* m}^{T}),

(89)

where

K_{* m}

is the covariance matrix between a test configuration

R_{*}

and α-th inducing set

z_{α}

⁠. Furthermore, we weight the prediction with

β_{α}

following the concepts of the robust BCM where the individual committee prediction is weighted by the differential entropy,

β_{α} = - log (σ_{α}^{2})

⁠. In the BCM, a particular local expert, which has the highest similarity with a test configuration

R_{*}

⁠, mostly contributes to the energy, forces, and pressure predictions. The final normalization factor

\hat{S}

is thus given as

\hat{S} = {[\sum_{α}^{P} \frac{β_{α}}{σ_{α}^{2}}]}^{- 1} .

(90)

The key advantage of using the SGPR-based BCM is that it not only provides sparsification but also allows us to invert an approximate $(m / P)$ -dimensional matrix rather than a full m-dimensional one (⁠ $m = \sum_{α}^{P} m_{α}$ ⁠). This reduces the computational cost of the matrix inversion from $O (m^{3})$ to $O (m^{3} / P^{3})$ ⁠. Compared to “product of experts” (POE), where the predicted potential energy is given by $E \approx \frac{1}{P} \sum_{α}^{P} E_{α}$ ⁠, the BCM employs a weighting factor $β_{α} / σ_{α}^{2}$ to ensure that models with higher confidence contribute more significantly to the predicted energy compared to less confident models.

It is validated that the BCM potential accurately estimates the melting temperatures of ice II, III, and V under pressure conditions of 0.2–0.6 GPa.¹¹² To achieve this, on-the-fly active learning SGPR NpT MD simulations was performed to create BCM MLP for each coexisting system: ice II-water, ice III-water, and ice V-water. Three separately trained SGPR local expert MLP models were combined into one sparse robust BCM potential. This combined potential was used for NpH MD simulations of coexisting ice–liquid systems. Non-active NpH MD simulations were performed to determine the melting temperatures of ice II, III, and V. During these MD simulations, the conservation of enthalpy (H) was ensured while monitoring fluctuations in total energy, which includes both potential and kinetic energy.

Figure 8 illustrates the melting points of ice II, III, and V as determined from our NpH MD simulations using the MB-pol water model. The ice–water phase diagram predicted with the MB-pol potential closely aligns with those obtained from enhanced-coexistence simulations, as reported by Bore and Paesani.¹¹³ The BCM MLP stabilizes ice V over ice II, consistent with Zhang et al.'s predicted phase diagram.¹¹⁴ However, using a sophisticated deep NN potential, Bore and Paesani¹¹³ observed that ice V became a metastable phase, absent from the phase diagram. This discrepancy is attributed to the failure in incorporating high model uncertainty configurations of the water–ice system in the training dataset. In contrast, the active BCM potential adeptly learns model uncertainties on-the-fly, enhancing model accuracy and thereby stabilize ice V in the phase diagram (Fig. 8). However, the melting points predicted by the BCM MLP exhibit a slight upward shift in temperature relative to the benchmark MB-pol potential (Fig. 8). This deviation primarily stems from the noise present in the energy and force predictions of the BCM potential. Notwithstanding the shifts in melting temperatures, we emphasize the BCM potential's capability to accurately replicate the physical and thermodynamic properties of the reference MB-pol water model. While the discrepancies in melting temperatures suggest a room for further improvements in accuracy of model, the significance lies in the BCM potential's reproduction of reference phase diagram. This affirms the BCM potential's utility in capturing the essence of water's behavior under various conditions while reducing the computational scaling of training SGPR MLP for large data.

FIG. 8.

View large Download slide

Ice–water phase diagram. Ice–water melting points are calculated through NpH MD simulations of coexisting ice–liquid systems using the MB-pol (circle lines) and BCM (star dashed lines). Solid black line represents the experimental melting line.¹¹⁷ Figure taken from Ref. 112. Reproduced with permission from Willow et al., Phys. Chem. Chem. Phys. 26, 22073–22082 (2024). Copyright 2024 Royal Society of Chemistry.

Q. Comparison of GNNs and SGPR

Figure 9 compares GNNs and SGPR during training and testing. The Material Project Trajectory (MPTrj)⁹⁷ serves as an example of a large dataset for training, containing about 1 500 000 compounds. The NNs^95–98 used the entire dataset directly, employing an embedding layer, interaction layers, and a readout layer to predict output values such as energy, with optimization. In contrast, SGPR organizes the data into subsets based on the similarity of LCEs, creating the local expert models. Each local expert is built from its unique LCEs, which serve as an inducing dataset and a solution (w) to predict energy. These local expert models form the basis for built-in SGPR kernels. To estimate the potential energy of a test configuration $x^{*}$ (⁠ ${ρ_{i}^{*}}$ ⁠), NNs rely on the trained parameters, but SGPR selects local experts (⁠ $h, j, k)$ and quantifies uncertainty based on the similarity measure $k (ρ_{i}^{*}, ρ_{j}^{g})$ between the LCE (⁠ $ρ_{i}^{*}$ ⁠) of the testing configuration and the existing inducing LCEs ${ρ_{j}^{g}}_{j = 1}^{m}$ ⁠.

FIG. 9.

View large Download slide

Comparison between GNNs and SGPR during training and testing. Multiple stacked layers symbolize a large dataset covering materials such as ice/water, batteries, perovskites, zeolites, metal–organic frameworks (MOFs), molten salts, and heterogeneous catalysis. Each stack corresponds to configurations with the same atomic numbers. The GNN model consists of an embedding layer, interaction layers, and a readout layer for predicting output values. During training, node features for all atomic numbers are initialized through the embedding layer and updated via interaction layers based on the edge information. The readout layer then converts the updated node features into the output values, in particular energy values. All model weights are optimized using the Adam optimizer. During testing, the train modeled model predicts the energy of a test configuration, labeled as $x^{*}$ in a green box. For SGPR, during training, configurations in the dataset are classified into subsets based on their similarity. Yellow boxes labeled “expert” refer to SGPR local expert model, which includes an inducing set (z) and a solution $(w)$ ⁠, as described in Eqs. (74) and (77). In testing, the predicted energy value $E^{*}$ is obtained from the combination of experts' kernels (⁠ $h, j, k$ ⁠) keeping the similarity $k (ρ_{i}^{*}, ρ_{j}^{g})$ between LCEs $ρ_{i}^{*}$ in $x^{*}$ and those $ρ_{j}^{g}$ in inducing sets z. For one universal expert model by combining these local experts, we use POE, MOE, or FOE as described in Sec. II O as well as BCM as described in Sec. II P.

For uncertainty quantification (UQ), SGPR estimates uncertainty based on the similarity measurement $k (ρ_{i}^{*}, ρ_{j}^{g})$ ⁠, providing UQ values [Eqs. (79) and (89)] using similarity measures, particularly in regions far from the training data, between the LCE of the testing configuration and the existing inducing LCEs. NNs, on the other hand, face difficulties in producing robust UQ estimates.^115,116 Methods such as Bayesian NNs, including Bayes by Backprop,^118,119 Monte Carlo-Dropout,^120,121 and deep ensembles,¹²² have been proposed to estimate uncertainties in the predicted total energy, but not for local energy predictions. Other UQ approaches use feature space,¹²³ latent space,¹²⁴ or GP-based kernel on latent space,¹²⁵ with the Euclidean norm distance in latent features showing some success. Recently, UQ could be improved by analyzing variations in predicted forces rather than predicted energies,¹²⁶ while handling variations in predicted forces remain less reliable in data-sparse regions. In this regard, SGPR stands out for its principled approach, providing closed-form uncertainty estimates, data-dependent variance, and the flexibility of a non-parametric model, which adapts its complexity to the data.

As for scalability, GNNs are able to handle large datasets and many atoms efficiently through local message-passing mechanisms. However, significant disadvantages of GNN models arise from the significant training data requirements, UQ, interpretability, and computational complexity during training. GNNs typically require large datasets to reach optimal performance. If training data are limited, they may not perform as well as GPR-based models, which are better at handling small data regimes. The need for extensive data can be a drawback in situations where high-quality labeled data (e.g., quantum mechanical calculations) are expensive or time-consuming to generate. GNNs are generally deterministic models, without providing UQ. This makes them less useful for tasks where knowing the confidence in prediction is crucial. GNNs are often seen as black-box models, making it difficult to interpret why certain predictions are made. This lack of transparency can be a concern in certain applications where interpretability is important. The underlying physics captured by the model can be harder to analyze compared to GPR where the kernels can often be mapped to physical processes. Although GNNs are scalable during inference, training them can still be computationally expensive, especially for equivariant models that require careful enforcement of rotational symmetries and use large datasets.

GPR models typically rely on SOAP descriptors, which become computationally intensive as the number of elements increases. This is because the descriptor needs to account for the diverse chemical environments and interactions for each additional element. The computational cost in GPR grows significantly when a greater number of atomic species is involved, making it less scalable for systems with many different elements. On the other hand, GNN models use graph-based embeddings, which capture structural information about the molecules without scaling as severely with the number of elements. GNNs are more flexible and efficient for large systems due to the ability to generalize molecular structures through learned embeddings. GNNs are highly scalable by handling large datasets and systems with many atoms efficiently due to their local message-passing nature. However, significant disadvantages of GNN models arise from the significant training data requirements, UQ, interpretability, and computational complexity during training. GNNs typically require large datasets to reach optimal performance. If training data are limited, they may not perform as well as GPR-based models, which are better at handling small data regimes. The need for extensive data can be a drawback in situations where high-quality labeled data (e.g., quantum mechanical calculations) are expensive or time-consuming to generate. GNNs are generally deterministic models, without providing UQ. This makes them less useful for tasks where knowing the confidence in prediction is crucial. GNNs are often seen as black-box models, making it difficult to interpret why certain predictions are made. This lack of transparency can be a concern in certain applications where interpretability is important. The underlying physics captured by the model can be harder to analyze compared to GPR where the kernels can often be mapped to physical processes. Although GNNs are scalable during inference, training them can still be computationally expensive, especially for equivariant models like NequIP and MACE that require careful enforcement of rotational symmetries and use large datasets.

On the other hand, GPR is a non-parametric Bayesian approach used for regression tasks. SGPR is a variant of GPR that reduces the computational cost by selecting a subset as a natural way to quantify uncertainty in the predictions. This is useful in materials discovery or simulations where understanding the reliability of predictions is critical. The Bayesian nature of SGPR allows for the prediction of both mean and variance, which is a unique strength over many GNNs. SGPR models tend to perform well with relatively small and moderate datasets, as they rely on a kernel-based approach that implicitly embeds prior knowledge about smoothness in the data. SGPR is particularly useful when training data are scarce or expensive to obtain. GPR models are highly interpretable, as the choice of kernels can be directly related to physical principles, and the model predictions can be understood as weighted averages of the training data. The ability to choose specific kernel functions provides flexibility to capture particular characteristics of atomic interactions. The predictions of GPR are generally smooth and conservative, which can be beneficial in certain applications, especially when physical smoothness constraints are essential, like for energy landscapes. In contrast, there are disadvantages of GPR. As the number of atoms and interactions increases, even sparse approximations become computationally extensive for very large systems. GPR relies heavily on the choice of kernel function to model interactions. To mitigate the scaling issue of GPR, feature compression, sparse representations, and partitioning strategy could be exploited. Techniques such as dimensionality reduction on SOAP descriptors (e.g., principal component analysis) could help reduce the complexity of the feature space without losing essential chemical information. Using a sparser representation of the SOAP descriptor could lessen the computational load by focusing on the most relevant atomic environments. Breaking down the chemical space into smaller regions (e.g., clustering molecules based on their element types) may allow the model to focus on subsets of elements at a time.

Overall, GNNs excel in scalability, can model complex interactions, and are better suited for large datasets, but require more data and computational resources, and lack built-in uncertainty quantification. SGPR is an excellent choice for small datasets, interpretable models, and applications where uncertainty quantification is critical. While GPR struggles with scalability and expressiveness in complex systems, SGPR mitigates such issues and, furthermore, the built-in-library SGPR potential/force-fields would maintain the accuracy of complex materials model systems without significant loss of the computational speed, enhance the prediction power of structure, energetics, and properties of unknown novel materials, and help in understanding of their physical science concepts.

There is no definitive answer for the exact maximum number of elements SGPR can handle, as it depends on the specific implementation and computational resources available. However, most SGPR models may struggle beyond 6–10 distinct elements, especially when high accuracy and rich descriptors are required. Meanwhile, GNNs are known for their strong extrapolation capabilities due to the generality of their learned embeddings. In contrast, GPR-based models are generally more interpolative, meaning they excel at making predictions within the bounds of the data they were trained on. While GPR models show limitations in scalability and extrapolation, particularly when the number of elements grows, GNN models provide a more robust framework for handling diverse chemical systems. In particular, SGPR-based models offer significant advantages in smaller, moderate-sized well-defined systems where interpretability and fine-grained control over descriptors like SOAP are essential. Overall, the choice between SGPR and GNN-based methods depends on specific application needs, data size, and computational resources. SGPR is preferred for smaller datasets requiring interpretable models and theoretically grounded uncertainty estimates, while GNNs are favored for large datasets and scalable predictions. SGPR provides a more interpretable and theoretically grounded uncertainty estimate, while NN-based methods are scalable but often rely on approximations or ensembles that may be less reliable in data-sparse regions. In this regard, scalable built-in-library SGPR could be a useful approach for moderately larger datasets. Future research may explore ways to combine the strengths of both SGPR and GNN approaches for broader applicability, particularly for larger datasets in materials science.

R. Summary of Sec. II

In Sec. II, we discussed ML representations of the ab initio PES. Two of the most important branches are NN³⁶ and GAPs.⁴⁷ NNs are computationally efficient because they scale as $O (n)$ with the size of the data. However, since they contain many hyper-parameters, they are more likely to over-fit the training data and not generalize well. Therefore, NNs are often used when large datasets are available. On the other hand, Gaussian process based methods are almost non-parametric; only a few hyper-parameters are present in the covariance kernel. They work well even with small datasets, but their computational cost grows fast with the size of data. The cost for building the covariance matrix scales as $O (n^{2})$ and for its inversion as $O (n^{3})$ ⁠. Therefore, approximate inference methods for the GPR, which transform these unfavorable cost scaling to $O (n)$ ⁠, are crucial. BLR method^43,91) is one of the possible solutions. However, a method better than linear regression is sparse approximation.

Another important issue is on-the-fly modeling and active learning^43,44,91,101 where the data required for building the model is generated alongside molecular dynamics (MD) simulations. For this purpose, Gaussian process based regression methods are much more suitable than NN because the predicted variance can be used as a proxy for on-the-fly accuracy of the model and a criterion for data generation by ab initio calculations. Although the uncertainty can be emulated with NN by training identical networks with different initial random seeds, updating an existing Gaussian process model is much easier.

Currently, Gaussian process based active learning algorithms use the predictive variance [see Eq. (40)] as a criterion for on-the-fly data generation. A common objective for these methods is to show that their predictive variance correlates well with the actual ML errors. However, we noticed that the predictive variance is blind to the energy and forces data. For this, we devised a new criterion that includes these data as well. As a result, during on-the-fly sampling, the model is generated with much smaller data while maintaining or improving the accuracy, which dramatically increases the efficiency. This is accomplished by low-rank approximation of the full covariance matrix in normal GPR. Using SGPR, the computational cost is effectively reduced to that of the linear regression methods while maintaining the accuracy. In addition, we discussed a novel adaptive sampling algorithm. The excellent performance of this formalism (SGPR + adaptive sampling) is demonstrated with benchmark simulations. As a benchmark simulation, we applied this algorithm to molecular dynamics simulations of bulk Si. The SGPR potential, which was generated using only 13 ab initio calculations, regenerated the phonon spectra almost identical to DFT. The combinatorial approach to generate a more universal expert system dealing with more elements from several expert systems dealing with less elements opens the door for strategies like “high level parallelism” or “a library of SGPR potentials,” which can be combined by simple cross covariance calculations (based on the MOE approach) for the cases to describe increasingly complex materials. Therefore, this approach can be useful for modeling wide-ranging composite materials for tailoring the performance. Since the computational cost of SGPR is proportional to the number of atoms N in the system, much bigger systems can be efficiently simulated compared with DFT where the computational cost scales as $\sim N^{3}$ ⁠. This allows us to study large and complex systems.

III. BASIC FUNCTIONS

The kernel regression methods in Sec. II assumed that a latent energy function $E (d_{i})$ exists for the contribution of atoms i in the potential energy of the system, where $d_{i} = D (ρ_{i})$ is the descriptor for the local chemical environment of atom i. Then using an appropriate covariance kernel $K (d_{i}, d_{j}) = 〈 E (d_{i}) E (d_{j}) 〉$ ⁠, representations based on linear regression or Gaussian process regression are found for $E$ using a set of ab initio data. Here we discuss the kernels and descriptors in more details.

A. General aspects

In all of the previous derivations, we assumed that the latent energy function (LEF) is related to the potential energy by Eq. (3), where every atom i has its own descriptor d_i. Although this can be applied everywhere, it is not the only possible approach. For instance, a group of atoms can be described with a single descriptor and the associated LEF. This is specially useful for molecular systems where coordinates of a group of atoms are strongly correlated. In such cases, the generalized expression in Eq. (1) can be assumed. The regression equations can still be applied with simple adjustments: ρ_i represents the translationally invariant degrees of freedom in the group of atoms indexed by i.

Ideally, the descriptor should be invariant with respect to the symmetry operations that do not change the potential energy. Symmetries that we are most concerned with are translations, rotations, and permutations of identical atoms. In some cases, certain symmetries can be ignored in the descriptor design and enforced manually. For instance, in a system with strong chemical bonds that permutations of atoms never occur, one can ignore this symmetry. Instead, the indices that are chosen for the atoms should be kept fixed.

In designing the descriptors for LCEs, since ρ_i defined by Eq. explicitly depends on the cutoff radius r_c, once must make sure that the descriptor function

D (ρ_{i})

is continuous and differentiable with respect to the coordinates of atoms in the vicinity of r_c

lim_{| r_{i j}^{μ} | \to r_{c}^{-}} \frac{\partial D (ρ_{i})}{\partial r_{i j}^{μ}} = lim_{| r_{i j}^{μ} | \to r_{c}^{+}} \frac{\partial D (ρ_{i})}{\partial r_{i j}^{μ}} = 0 .

(91)

The above behavior can be enforced by appropriate inclusion of smooth cutoff functions in the descriptor. The simplest choice for the cutoff function is

Θ (| {\vec{r}}_{i j} |) = {\begin{matrix} {(1 - \frac{| {\vec{r}}_{i j} |}{r_{c}})}^{2} & | {\vec{r}}_{i j} | \leq r_{c}, \\ 0 & | {\vec{r}}_{i j} | > r_{c} . \end{matrix}

(92)

Another commonly used cutoff function is³⁶

Θ (| {\vec{r}}_{i j} |) = {\begin{matrix} \frac{1}{2} (cos (\frac{| {\vec{r}}_{i j} |}{r_{c}} π) + 1) & | {\vec{r}}_{i j} | \leq r_{c}, \\ 0 & | {\vec{r}}_{i j} | > r_{c} . \end{matrix}

(93)

B. Combining descriptors

The LEF can be split into several parts, for instance, to two-body (2b), three-body (3b), and many-body (mb), where each has its own descriptor. Formally,

d_{i} = \underset{α}{\otimes} d_{i}^{α},

(94)

E (d_{i}) = \sum_{α} E^{α} (d_{i}^{α}),

(95)

where α represents the index of different types of descriptors. Then for the covariance kernel, we can write

K (d_{i}, d_{j}) = \sum_{α β} 〈 E^{α} (d_{i}^{α}) E^{β} (d_{j}^{β}) 〉 .

(96)

With the assumption that

〈 E^{α} (d_{i}^{α}) E^{β} (d_{j}^{β}) 〉

is zero unless

α = β

⁠, the kernels become additive

K^{α} (d_{i}^{α}, d_{j}^{α}) = 〈 E^{α} (d_{i}^{α}) E^{α} (d_{j}^{α}) 〉,

(97)

K (d_{i}, d_{j}) = \sum_{α} K^{α} (d_{i}^{α}, d_{j}^{α}) .

(98)

C. Rotational symmetry

In describing the LCEs, we are interest in quantities that remain invariant with respect to rotations. Any function of the type

f (θ, ϕ)

⁠, where θ and

ϕ

are polar and azimuthal angles in spherical coordinates, respectively, can be represented as

f (θ, ϕ) = \sum_{l = 0}^{\infty} \sum_{m = - l}^{l} c_{l m} Y_{l}^{m} (θ, ϕ),

(99)

where

Y_{l}^{m}

are the spherical harmonics. Therefore, it suffices to know the transformations of the spherical harmonics under applied rotations. If

(θ, ϕ) \to (θ', ϕ')

under rotation R, then

Y_{l}^{m} (θ', ϕ') = \sum_{m' = - l}^{l} {[D_{m m'}^{l} (R)]}^{*} Y_{l}^{m'} (θ, ϕ),

(100)

where

D_{m' m}^{l} (R)

are elements of the Wigner D-matrix, which is a block-diagonal matrix with the property

D_{m' m}^{l} (R^{- 1}) = {[D_{m m'}^{l} (R)]}^{*} .

Equation in matrix notations (⁠

Y_{l}

is a column vector of the length

2 l + 1

⁠, and

D^{l}

a square matrix of the shape

(2 l + 1) \times (2 l + 1)

⁠) becomes

Y_{l} (θ', ϕ') = {[D^{l} (R)]}^{*} Y_{l} (θ, ϕ) .

(101)

Considering that

{(D^{l})}^{†} D^{l} = I^{l}

⁠, it is easy to see that

Y_{l}^{†} (θ, ϕ) Y_{l} (θ, ϕ)

remains invariant with respect to rotations. Thus if we define

f_{l} = {[c_{l, - l} Y_{l}^{- l} (θ, ϕ), \dots, c_{l, l} Y_{l}^{l} (θ, ϕ)]}^{T}

⁠, q_l defined below is invariant under rotations such that

q_{l} (θ, ϕ) = \sum_{m = - l}^{l} {[c_{l m} Y_{l}^{m} (θ, ϕ)]}^{*} c_{l m} Y_{l}^{m} (θ, ϕ) = f_{l}^{†} f_{l} .

(102)

In fact, this is used in defining the bond orientational order parameter, which has been used as the descriptor. More generally,

Y_{l}^{†} (θ_{1}, ϕ_{1}) Y_{l} (θ_{2}, ϕ_{2})

is invariant with respect to rotations. The following relations will be useful:

\int d R D_{μ ν}^{λ} (R) {[D_{m m'}^{l} (R)]}^{*} = \frac{8 π^{2}}{2 l + 1} δ_{l λ} δ_{m μ} δ_{m' ν},

(103)

\int d R D_{μ ν}^{λ} (R) D_{β γ}^{α} (R) {[D_{m m'}^{l} (R)]}^{*} = \frac{8 π^{2}}{2 l + 1} C_{λ μ, α β}^{l m} C_{λ ν, α γ}^{l m'},

(104)

where

C_{λ μ, α β}^{l m}

are the Clebsch–Gordan coefficients and the integration over all rotations can be done using Euler angles, i.e.,

\int d R = \int_{0}^{2 π} d α \int_{0}^{2 π} d γ \int_{0}^{π} sin β d β

⁠.

D. Two-body and many-body

Application of two-body (2b) descriptor is straightforward since distances are invariant with respect to rotations or permutations of atoms. If the RBF kernel is chosen as the covariance kernel, then 2b kernel becomes

k^{2 b} (ρ_{i}, ρ_{j}) = \sum_{k \in ρ_{i}, l \in ρ_{j}} exp (- \frac{{(| {\vec{r}}_{i k} | - | {\vec{r}}_{j l} |)}^{2}}{2 ξ^{2}}) Θ (| {\vec{r}}_{i k} |) Θ (| {\vec{r}}_{j l} |),

(105)

where Θ is a smooth cutoff function as Eq. or . For a three-body (3b) descriptor, one has to construct an array that is invariant with respect to the permutation of atoms. A possible choice is⁴⁸

d_{i}^{j k} = [| {\vec{r}}_{i j} | + | {\vec{r}}_{i k} |, {(| {\vec{r}}_{i j} | - | {\vec{r}}_{i k} |)}^{2}, | {\vec{r}}_{j k} |],

(106)

where atoms j and k are present in the environment of atom i. The above array is invariant with respect to

j \leftrightarrow k

permutation.

E. Coulomb matrix

Another descriptor that has been widely used is Coulomb matrix,¹²⁷ which is defined by

M_{i j} \propto \frac{Θ (| {\vec{r}}_{i j} |)}{| {\vec{r}}_{i j} |},

(107)

where Θ is a optional smooth cutoff function. M_ij is not invariant with respect to permutations of atoms. Therefore, a fingerprint is often generated using the sorted eigenvalues of the Coulomb matrix

d_{i} = [λ_{1}, λ_{2}, λ_{3}, \dots] λ_{1} \geq λ_{2} \geq λ_{3} \geq \dots .

(108)

Still, discontinuities may occur if a level crossing happens in the eigenvalues of the Coulomb matrix. This descriptor is often used for studies of molecular systems.^50,51

F. Behler–Parrinello symmetry functions

A many-body descriptor for the LCEs, which is often used in neural networks, is Behler–Parrinello symmetry functions.³⁶ Radial symmetry functions are defined by

g_{i}^{1} (η, r_{s}) = \sum_{j \in ρ_{i}} e^{- η {(| {\vec{r}}_{i j} | - r_{s})}^{2}} Θ (| {\vec{r}}_{i j} |)

(109)

and angular symmetry functions by

g_{i}^{2} (η, λ, ζ) = 2^{1 - ζ} \sum_{j k \in ρ_{i}} {(1 + λ cos θ_{ijk})}^{ζ} e^{- η (| {\vec{r}}_{i j} |^{2} + | {\vec{r}}_{i k} |^{2} + | {\vec{r}}_{j k} |^{2})} Θ (| {\vec{r}}_{i j} |) Θ (| {\vec{r}}_{i k} |) Θ (| {\vec{r}}_{j k} |),

(110)

where

θ_{ijk} = \frac{{\vec{r}}_{i j} . {\vec{r}}_{i k}}{| {\vec{r}}_{i j} | | {\vec{r}}_{i k} |}

(111)

and

ω = (η, r_{s}, λ, ζ)

are the free parameters. The symmetry functions with different values of ω are packed into a descriptor vector as

d_{i} = [(g_{i}^{1} (ω_{1}), g_{i}^{2} (ω_{1})), (g_{i}^{1} (ω_{2}), g_{i}^{2} (ω_{2})), \dots] .

(112)

G. Smooth overlap of atomic positions

For a similarity/covariance kernel between two LCEs, the SOAP is defined as¹²⁸

k (ρ_{i}, ρ_{i'}) = \int d R | \int d \vec{r} ϱ_{i} (\vec{r}) ϱ_{i'} (R \vec{r}) |^{2},

(113)

where

R

is the 3D rotations operator,

ϱ_{i} (\vec{r})

is the atomic density defined as

ϱ_{i} (\vec{r}) = \sum_{j \in ρ_{i}} exp (- α | \vec{r} - {\vec{r}}_{i j} |^{2}),

(114)

and α is a hyper-parameter. Such a kernel is invariant with respect to all rotations

k (R_{1} ρ_{i}, R_{2} ρ_{i'}) = k (ρ_{i}, ρ_{i'}) .

(115)

The above integrals can be carried out analytically, which ignoring a constant multiplier, results in

k (ρ_{i}, ρ_{i'}) \propto \sum_{l, m, m'} \frac{1}{2 l + 1} I_{m m'}^{l}^{*} I_{m m'}^{l},

(116)

where

I_{m m'}^{l} = \sum_{j \in ρ_{i}, j' \in ρ_{i'}} e^{- α (r_{j}^{2} + r_{j'}^{2}) / 2} ι_{l} (α r_{j} r_{j'}) Y_{l}^{m} ({\hat{n}}_{j}) Y_{l}^{m'}^{*} ({\hat{n}}_{j'}),

(117)

r_{j} = | {\vec{r}}_{i j} |, {\hat{n}}_{j} = {\vec{r}}_{i j} / r_{j}

⁠, and ι_l is the modified Bessel function of the first kind (Bessel-i). Unfortunately, these equations are not computationally efficient because of the dependence of

I_{m m'}^{l}

identities on both environments. In the study by Bartók et al.,¹²⁸ this is remedied by expansion of the atomic density using an orthonormal radial basis set

{g_{n}}

ρ (\vec{r}) = \sum_{n, l, m} c_{nlm} g_{n} (r) Y_{l}^{m} (\hat{n}),

(118)

which transforms the

I_{m m'}^{l}

identities into inner products

I_{m m'}^{l} = \sum_{n} c_{nlm} c_{nlm'}^{' *},

(119)

which in turn transforms the kernel into an inner-product

k (ρ, ρ') \propto \sum_{n, n', l} (\sum_{m} \frac{c_{nlm} c_{n' l m}^{*}}{\sqrt{2 l + 1}}) (\sum_{m'} \frac{c_{nlm'}^{' *} c'_{n' l m'}}{\sqrt{2 l + 1}}) = \sum_{n, n', l} p_{n n' l} p_{n n' l}^{' *} .

(120)

Here we achieve the same goal, by considering the series expansion of Bessel-i functions ι_l,¹²⁹

ι_{l} (α r_{j} r_{j'}) = \sum_{n = 0}^{\infty} a_{n l} r_{j}^{2 n + l} r_{j'}^{2 n + l},

(121)

where

a_{n l} = {(α / 2)}^{2 n + l} / (n + l)! n!

⁠. Using this expansion, we get the same expression for

I_{m m'}^{l}

as Eq. with the expansion coefficients already determined:

c_{nlm} = \sqrt{a_{n l}} \sum_{j \in ρ_{i}} e^{- α r_{j}^{2} / 2} r_{j}^{2 n + l} Y_{l}^{m} ({\hat{n}}_{j}) .

(122)

For elimination of discontinuities due to the application of a cutoff, we need to modify the descriptor, for instance, by the substitution

e^{- α r_{j}^{2} / 2} \to Θ (r_{j}) e^{- α r_{j}^{2} / 2}

⁠, where Θ is a smooth cutoff function with a continuous derivative at r_c. The a_nl multipliers decrease extremely fast by increasing n and l; therefore, it suffices to calculate up to only a few orders of expansion. α = 1, r_c = 8, and

l, n \leq 3

are chosen for reasonable accuracy.

Extension to multi-species environments is straightforward.¹³⁰ If the coefficients c_nlm are calculated for each atomic species in the LCE, the total descriptor vector is

p_{n n' l}^{α β} = \frac{1}{\sqrt{2 l + 1}} \sum_{m} c_{nlm}^{α} c_{n' l m}^{β}^{*},

(123)

where α and β denote the species. Then, the kernel may be generalized to

k (ρ_{i}, ρ_{i'}) = μ_{s_{i} s_{i'}} {(\sum_{α, α', β, β^{'}, n, n', l} {[p_{n n' l}^{α β}]}_{i} {[p_{n n' l}^{α' β'}]}_{i'} χ_{α α'} χ_{β β'})}^{η},

(124)

where

χ_{α β}

is an alchemical similarity kernel between the atomic species in LCEs and

μ_{s_{i} s_{i'}}

is a kernel that depends on the species of the central atoms

i, i'

and quantifies the covariance between local energies. Here we have chosen

χ_{α β} = δ_{α β}

and

μ_{s_{i} s_{i'}} = δ_{s_{i} s_{i'}}

⁠. The latter corresponds to summation of independent Gaussian processes defined for each atomic species. The exponent η controls the sensitivity of the kernel, which is fixed at 4 in our work. Finally, the kernel can be normalized by

k (ρ_{i}, ρ_{i'}) \to \frac{k (ρ_{i}, ρ_{i'})}{\sqrt{k (ρ_{i}, ρ_{i}) k (ρ_{i'}, ρ_{i'})}} .

(125)

H. Summary of Sec. III

In this section, we discussed combinations of covariance kernels such as two-body, three-body, and many-body. As examples of many-body descriptors, we represented Coulomb matrix, Behler–Parrinello symmetry functions, and smooth overlap of atomic positions (SOAP). We found analytical expressions for the components of the SOAP descriptor using the series expansion of the Bessel-i functions, which can be used for fast numerical evaluation of the descriptors. The SOAP kernel is almost non-parametric. Only a few parameters are present in the kernel, which have simple intuitive interpretations. Therefore, we can often skip iterative refinement using the gradient descent algorithm, which is computationally an exhaustive task. The SOAP descriptor with the series expansion of the Bessel-i functions (SE-SOAP) is used in this study.

IV. DYNAMIC PROPERTIES

In this section, we discuss diffusivity in solids and the methods used to estimate ionic conductivity.

A. Conductivity

Resistance R of a conductor with the length l and cross-sectional area A is given by

R = ρ \frac{l}{A},

(126)

where resistivity (or specific resistance) ρ has SI unit of Ωm. Conductivity (or specific conductance) is the reciprocal of resistivity

σ = \frac{1}{ρ} .

(127)

SI unit of conductivity is S/m, where Siemens S is 1/Ω.

B. Ionic conductivity in crystals

Solid electrolytes that allow fast ion conduction may lead to smaller and safer portable energy storage devices than the current batteries that come with liquid electrolytes.¹³¹ Cations such as Li⁺ diffuse in solid in the framework of anions or poly-anionic moieties. Many fast ion conductor solids have a body-centered-cubic arrangement of anions.¹³² In crystals, ionic conductivity rises from the formation and hopping of defects such as vacancies and interstitials whose concentration is temperature dependent. If defects are uncorrelated, conductivity can be described by an Arrhenius-type relationship

σ = qnu = σ_{0} T^{m} e^{- E_{a} / k_{B} T},

(128)

where

q, n, u, E_{a}

are the charge, concentration, mobility, and activation energy of mobile ions, respectively, and m is typically −1. E_a entails the energies required for the creation and hopping of mobile charges. In the simplest case where an ion jumps from one lattice site to a neighboring vacancy, we can write

σ_{0} = z \frac{n q^{2}}{k_{B}} e^{Δ S_{m} / k_{B}} α_{0}^{2} ν_{0},

(129)

where

z, Δ S_{m}, α_{0}, ν_{0}

are the geometric factor, entropy of migration, jumping distance, and attempt frequency, respectively. ν₀ is related to the character of lattice vibrations.

Diffusion is proportional to conductivity through the Nernst–Einstein relationship

D = \frac{u}{q} k_{B} T = \frac{σ}{n q^{2}} H_{R} k_{B} T,

(130)

where the Haven ratio H_R is the correlation factor that depends on the diffusion mechanism.¹³³ For uncorrelated ion hopping H_R = 1 and for highly correlated ionic motion observed in fast-ion conductors

H_{R} > 2

⁠. It has been suggested that concerted migrations of multiple ions with low energy barriers is the origin of fast ionic diffusion in superionic conductors.¹³⁴

C. Diffusion in solids

Nudged elastic band (NEB) calculations are used to calculate the energy barrier along a specific diffusion path if dominant diffusion pathways are known. For a concerted diffusion mechanism, this becomes even more complicated as the number of possible pathways become increasingly large. The ab initio molecular dynamics (AIMD) simulations, on the other hand, are able to realistically determine the diffusion coefficient and the conductivity of materials as well as the diffusion pathways. However, due to the huge computational resources required for the electronic structure calculations, these simulations are often applied only to small systems (⁠ $\sim 100$ atoms) and short intervals of time (⁠ $\sim 100 ps$ ⁠). Therefore, sometimes classical MD simulations that use empirical potentials or force fields are utilized to scale up calculations.

Another important aspect of the studies of diffusion in the solids is that, due to the Arrhenius-type temperature dependence of the diffusion events, they occur too slowly at room temperature to be sampled with statistical certainty using AIMD¹³⁵ even in the fastest ion conductors. Therefore, often simulations are carried out at elevated temperatures and the diffusion coefficient at room temperature is approximated by extrapolation. However, this is not possible in the cases where the system of interest is unstable at higher temperatures. In fact, such a case is encountered in the current study.

D. Linear response theory

The ionic conductivity σ is given by the Green–Kubo relation in linear response theory¹³⁶

σ = \frac{q^{2}}{V k_{B} T d} \int_{0}^{\infty} 〈 \vec{J} (t) \cdot \vec{J} (0) 〉 d t,

(131)

where d is the dimension, V is the volume, N is the number of mobile ions, q is the charge of ions, k_B is the Boltzmann constant, T is temperature, and

\vec{J} (t) = \sum_{i = 1}^{N} {\vec{v}}_{i} (t)

is the total particle flux. This can be written in the form of Nernst–Einstein equation as

σ = \frac{N q^{2}}{V k_{B} T} D_{σ},

(132)

where

D_{σ}

is called the jump diffusion coefficient and is related to the diffusion of the center of mass. The tracer diffusion coefficient is

D = \frac{1}{N d} \int_{0}^{\infty} \sum_{i = 1}^{N} 〈 {\vec{v}}_{i} (t) \cdot {\vec{v}}_{i} (0) 〉 d t .

(133)

One can see that

D_{σ} = D + \frac{1}{N d} \int_{0}^{\infty} \sum_{i \neq j} 〈 {\vec{v}}_{i} (t) \cdot {\vec{v}}_{j} (0) 〉 d t .

(134)

With the definition of Haven ratio as

H_{R} = D / D_{σ}

⁠,¹³⁷ ionic conductivity is often expressed in terms of the tracer diffusion coefficient as

σ = \frac{N q^{2}}{V k_{B} T} \frac{D}{H_{R}} .

(135)

If the motions of ions are uncorrelated, then H_R = 1.

E. Diffusion constant from molecular dynamics

The diffusion coefficient can be calculated from molecular dynamics simulations. The mean squared displacement (MSD) is defined as

MSD (t) = \frac{1}{N} \sum_{i = 1}^{N} | Δ {\vec{r}}_{i} (t) |^{2},

(136)

where N is the number of mobile ions and

Δ {\vec{r}}_{i} (t) = {\vec{r}}_{i} (t) - {\vec{r}}_{i} (0)

⁠. The tracer diffusion coefficient can be calculated from MSD using the relation

D = lim_{t \to \infty} \frac{1}{2 d t} 〈 MSD (t) 〉,

(137)

where d is the number of dimensions and

〈 〉

indicates the ensemble average. The squared mean displacement (SMD), which is related to the displacement of center of mass, is

SMD (t) = | \frac{1}{N} \sum_{i = 1}^{N} Δ {\vec{r}}_{i} (t) |^{2} .

(138)

Similarly, the jump diffusion coefficient can be obtained from SMD as

D_{σ} = lim_{t \to \infty} \frac{1}{2 d t} 〈 N SMD (t) 〉 .

(139)

Although, the diffusion of the center of mass is often too slow to be calculable with statistical certainty using AIMD.

F. Grain boundaries and long-range diffusivity

Consider an infinite 1D lattice where

p_{j} = exp (- ε_{j} / k_{B} T)

is the probability that a random walker will jump over the edge connecting j and j + 1 vertices (⁠

p_{j} = p_{j \to j + 1} = p_{j + 1 \to j}

⁠) and

ε_{j}

is the energy barrier. If p_j = p for all j, then

D = ν a^{2} p

is the diffusion coefficient for random walkers, where a is the lattice spacing and ν is the frequency of the jump trials. In general case where p_j is non-uniform, the diffusion coefficient becomes

D = ν a^{2} p^{f}

and

ε^{f} = - k_{B} T log (p^{f})

is the effective energy barrier, which can be calculated (analogous to series of resistors) from following:

\frac{1}{p^{f}} = 〈 \frac{1}{p_{j}} 〉 = 〈 exp (ε_{j} / k_{B} T) 〉 .

(140)

As a simple model for a crystal with grain boundaries, assume that one out of every l edge has a different energy barrier

ε_{g}

⁠; otherwise, the energy barrier is equal to

ε_{a}

⁠. The effective energy barrier becomes

ε^{f} = - k_{B} T log (\frac{l}{(l - 1) exp (ε_{a} / k_{B} T) + exp (ε_{g} / k_{B} T)}) .

(141)

If

\frac{ε_{a}}{k_{B} T}, \frac{ε_{g}}{k_{B} T}, \frac{1}{l} ≪ 1

⁠, then

ε^{f} \approx ε_{a} + \frac{ε_{g} - ε_{a}}{l},

(142)

which readily shows how the activation energy for long-range diffusion is related to the spacing of the grain boundaries.

In molecular dynamics simulations, often only the tracer diffusivity is calculated in a pure crystal. This can be used to estimate the short-ranged diffusivity. Grain boundaries and long-range diffusivity are difficult to consider because of the size limitations of ab initio simulations.

G. Diffusivity and activation energy

The diffusion coefficient (diffusivity) is defined as

D = lim_{t \to \infty} \frac{1}{2 d t} 〈 MSD (t) 〉,

(143)

where d = 3 is dimensions and

〈 \dots 〉

indicates the ensemble average. The directional MSD along

\hat{n}

is

{MSD}_{\hat{n}} (t) = \frac{1}{N} \sum_{i = 1}^{N} | Δ {\vec{r}}_{i} (t) \cdot \hat{n} |^{2} .

(144)

The activation energy E_a is calculated with the assumption of an Arrhenius temperature dependence for diffusivity as

D = D_{0} e^{- E_{a} / k_{B} T} .

(145)

The Nernst–Einstein relationship for conductivity σ is given by

σ = \frac{N q^{2}}{V k_{B} T} D,

(146)

where N is the number of the diffusing atoms, V is the volume, q is the ions electric charge, and T is the temperature. Alternatively, Arrhenius dependence for σT can be assumed and the diffusivity can be obtained by extrapolation. The numerical difference between the two cases is found to be small [within

O (1 %)

].

H. Summary of Sec. IV

Here we have described the dynamic properties such as ionic conductivity, diffusion, linear response theory, molecular dynamics simulations, grain boundaries, and long-range diffusivity for material simulations.

V. APPLICATIONS: SGPR-ASSISTED MATERIALS SIMULATIONS

A. Sulfide solid electrolytes

In this section, we discuss the Li diffusivity in sulfide solid electrolytes, which are one of the most promising candidates for future battery systems. Li₇P₃S₁₁ has shown high ionic conductivity. More powerful solid electrolytes showing remarkable ionic conductivity need to be developed for more promising battery systems. Here we investigate these materials using the SGPR algorithm and study its properties using large scale ML MD simulations. We also discuss the glass phase with a melt-quench simulation and the Li diffusivity in the glass phase and the crystallization.

FIG. 10.

Training data: for forces, RMSE is 0.13 eV/Å and R2=0.967; for energies, MAE is 3.817×10−3 eV per-atom.54 Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

View large Download slide

Training data: for forces, RMSE is 0.13 eV/Å and $R^{2} = 0.967$ ⁠; for energies, MAE is $3.817 \times 10^{- 3}$ eV per-atom.⁵⁴ Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

1. Li₂S-P₂S₅ system

We consider sulfide-based superionic conductors in the $(1 - z)$ Li₂S-zP₂S₅ system. The most stable compound in this system is considered to be Li₃PS₄. It exists as $γ$ -Li₃PS₄ below 195 °C and $β$ -Li₃PS₄ at higher temperatures. Depending on the experimental preparation method, the ionic conductivity as high as $\sim 10^{- 4} S / cm$ is reported for this material.¹³⁸ A much more conductive compound is the Li₇P₃S₁₁ glass-ceramic, which has a highly conductive phase due to the formation of a metastable crystal structure.¹³⁹ This crystal structure is found to be a triclinic cell made of PS₄ tetrahedra and P₂S₇ ditetrahedra.¹⁴⁰ With heat treatment to reduce the grain boundary resistance, the conductivity increases to $1.7 \times 10^{- 2} S / cm$ ⁠.¹⁴¹ Some of other fastest ion conductors are closely related to this system in that they are obtained upon substituting some of P or S atoms by other atoms (M, for example, Ge), namely, Li₂S-P₂S₅-GeS_n, Li₁₀GeP₂S₁₂,¹⁴² or a more complicated case of Li_9.54Si_1.74P_1.44S_11.7Cl_0.3.¹⁴³ Argyrodite type solid-state electrolytes (SSEs) are another family of electrolytes. Thus, the main degrees of freedom are to choose the ternary/quaternary systems, the composition ratios, and the route to the final product such as the heat treatment temperature and duration. First-principles theory calculations have been applied extensively for theoretical studies of solid electrolytes,¹⁴⁴ which helped in shedding light on the lithium's self-diffusion mechanism in Li₁₀GeP₂S₁₂,^145,146 the effects of anion/cation substitutions,¹⁴⁷ the cause of different ionic conductivity between tetragonal and cubic Li₇La₃Zr₂O₁₂,^148,149 and the effect of composition in Li₃OCl_1−xBr_x lithium-rich antiperovskite¹⁵⁰ and (1-z)Li₄SiO₄-(z)Li₃PO₄ solid solutions.¹⁵¹

2. Training

A sparse Gaussian process potential was trained with on-the-fly adaptive sampling algorithm for the well known crystal structures of γ-Li₃PS₄, β-Li₃PS₄, and Li₇P₃S₁₁.⁶³ For sampling data, the systems were quenched from ∼1500 to 700 K in ∼5 ps by velocity re-scaling and kept at this temperature for ∼45 ps using a time step of 2 fs. In quenching, there is no time for the crystal to melt (since atoms do not have time to displace too far the crystal becomes only locally disordered), but this local distortion results in stream of diverse LCEs to the sampling algorithm, which helps in faster stabilization of the model at the beginning of training. In total, only 456 ab initio calculations are performed from which only 111 are sampled by the model as the data X and 705 LCEs are sampled as the inducing set z.

3. Testing

The model on all 456 ab initio data generated during training fits ab initio forces with an RMSE of 0.13 eV/Å and a coefficient of determination of $R^{2} = 0.967$ (Fig. 10). To verify the model with independent data, ML MD simulations at five temperatures in the range 300–700 K for 0.8–2.2 ns are performed and 10³ snapshots are randomly selected. The accuracy of the model is tested with these samples in Fig. 11. For this testing set, we have RMSE = 0.14 eV/Å and $R^{2} = 0.944$ ⁠. The RMSE of the model is almost the same for the training data and testing sets. This shows that the model generalized very well.

FIG. 11.

Independent testing set: for forces, RMSE is 0.139 eV/Å and R2=0.944; for energies, MAE is 3.822×10−3 eV per-atom.54 Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

View large Download slide

Independent testing set: for forces, RMSE is 0.139 eV/Å and $R^{2} = 0.944$ ⁠; for energies, MAE is $3.822 \times 10^{- 3}$ eV per-atom.⁵⁴ Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

In addition, the liquid and glass phases of Li₇P₃S₁₁ are simulated with the SGPRs model. 50 + 70 samples are drawn from machine learning molecular dynamics (MLMD) simulation of the liquid phase Li₇P₃S₁₁ at 1200 K and glass phase (obtained after melt-quenching) at 600 K for ∼0.2 ns. The number of atoms in these samples is 672.

To show that the model is energy-conserving, a MD simulation was carried out in the NVE ensemble. The total energy is conserved up to small fluctuations of $\sim 10^{- 5}$ eV per atom, which are related to the finite time step of 2 fs. Since the model is trained at high temperatures, it is intriguing to see how well the model describes the properties close to 0 K. Specifically, NEB simulation is used in many ab initio studies to estimate the activation energy for Li diffusion since lengthy AIMD simulation is computationally very time consuming. We carried out NEB calculations with both the ML potential and direct DFT calculations. The trained potential in Fig. 12 describes NEB with an acceptable accuracy.

FIG. 12.

Nudged elastic band calculations for Li hopping to a neighboring vacancy (created by Li removal) in γ-Li3PS4 with DFT and ML potentials. (2×2×2) k-point grid is applied in DFT calculations.54 Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

View large Download slide

Nudged elastic band calculations for Li hopping to a neighboring vacancy (created by Li removal) in γ-Li₃PS₄ with DFT and ML potentials. $(2 \times 2 \times 2)$ k-point grid is applied in DFT calculations.⁵⁴ Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

4. Computational performance

For a system with N atoms, the computational cost of predicting energy and forces with the SGPR model is $O (N)$ ⁠, while the complexity of ab initio DFT calculations is $O (N^{3})$ or at best $O (N^{2} log N)$ unless a fragmental approach is not employed. The elapsed time for ab initio calculations of the energy and forces for a system with 672 atoms was 9718 s using VASP with 80 cores, while the elapsed time for the same calculation with SGPR was 5.3 s using only 20 cores. The cost defined by $elapsed \times cores / N$ is $\sim 0.15 s$ for SGPR, which is better than ab initio calculations by a factor of $O (10^{4})$ ⁠. This gain will grow exponentially as the number of atoms increases.

5. α and β phases

Large scale ML MD simulations of Li₇P₃S₁₁ (⁠ $4 \times 2 \times 2$ supercell, 672 atoms) were performed in isothermal-isobaric (NPT) ensemble at several temperatures in the range 300–1200 K and external pressure of $10^{5} Pa$ ⁠. For this size, the implementation of ML MD is $O (10^{4})$ faster than DFT MD. A phase transition was detected by ML MD (at $T \geq 450$ ⁠), which occurred by rotations of the P₂S₇ double-tetrahedra into a new orientational order. The initial and the new structures are referred to as α and β phases. The ground state energies of α and β phases obtained from exact DFT calculations are −4.413 and −4.416 eV per-atom, respectively, which shows that they are almost iso-energetic and the β phase is not an artifact of the SGPR potential. Also, examples from the β phase are encountered during training and testing. 0.14 eV/Å is the lumped RMSE for both phases. There appeared two steps in the MSD of S atoms corresponding to the delayed rotations of P₂S₇ double-tetrahedra in different layers. The delay was ∼0.5 ns at 450 K with a $4 \times 2 \times 2$ supercell. Specifically in the range 450–600 K (Fig. 13), the system remained a few hundred picoseconds in the α-phase before transition to the β-phase, which allows estimation of the Li diffusivity in both phases. In agreement with experimental reports,^139,140 at 900 K, the P and S atoms also started to diffuse, which indicates melting or its decomposition to more stable components.

From the diffusion coefficient (diffusivity) (Fig. 14) with the assumption of an Arrhenius temperature dependence for diffusivity, the activation energies E_a in α and β phases were 0.20 and 0.30 eV, respectively. The ionic conductivity in the α phase in NPT at RT(=300 K) becomes $3.5 \times 10^{- 2} S / cm$ ⁠, which is in reasonable agreement with previous DFT simulation results of $5.7 \times 10^{- 2}$ ¹⁵² in NVT and $4.5 \times 10^{- 2} S / cm$ ¹⁵³ within a limited box size. For the β phase, the RT conductivity is obtained by extrapolation of D and approximating V (which is T-dependent) by the volume of the relaxed cell. This yields the conductivity of $2.3 \times 10^{- 3} S / cm$ ⁠, which is more than an order of magnitude lower than the α phase. Alternatively, the Arrhenius dependence for σT is assumed. The numerical difference is $O (1 %)$ since the volume of α phase expands less than 1% (300–600 K) and β phase about $\sim 4 %$ (450–800 K).

The diffusivity was the largest along the $\vec{a}$ -axis in both phases. The main difference is that the diffusivity along $\vec{b}$ and $\vec{c}$ axes are almost equal in the α phase, while in the β phase, this symmetry is clearly broken. In Ref. 152, it is pointed out that the one-dimensional (1D) MSD of Li atoms along one of the axes in Li₇P₃S₁₁ is nearly a factor of 2 higher than that along the other two axes. This is consistent with the directional MSD.

6. Glass phase

The glass Li₇P₃S₁₁ is created by a melt-quench simulation from 1200 to 300 K. The Li conductivity of $9.3 \times 10^{- 3} S / cm$ is then obtained for the glass phase (with ML MD for 3 ns), which is higher than the conductivity in the β phase. Equilibration of the glass at higher temperatures shows that P atoms start to diffuse at 600 K, which indicates crystallization below this temperature (Fig. 13), in agreement with experiment.^139–141 The liquid and glass phases show RMSE = 0.21 eV/Å and $R^{2} = 0.917$ (see Fig. 15).

FIG. 13.

MSD of P atoms in glass phase Li7P3S11. Diffusion of these atoms at 600 K allows the glass to slowly relax and eventually crystallize. Figure taken from Ref. 54. Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

View large Download slide

MSD of P atoms in glass phase Li₇P₃S₁₁. Diffusion of these atoms at 600 K allows the glass to slowly relax and eventually crystallize. Figure taken from Ref. 54. Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

FIG. 14.

Li diffusivity (D) in α (blue) and β (red) phases of Li7P3S11. Figure taken from Ref. 54. Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

View large Download slide

Li diffusivity (D) in α (blue) and β (red) phases of Li₇P₃S₁₁. Figure taken from Ref. 54. Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

FIG. 15.

Liquid and glass phases: for forces, RMSE is 0.21 eV/Å and R2=0.917; for energies, MAE is 6.35×10−3 eV per-atom.54 Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

View large Download slide

Liquid and glass phases: for forces, RMSE is 0.21 eV/Å and $R^{2} = 0.917$ ⁠; for energies, MAE is $6.35 \times 10^{- 3}$ eV per-atom.⁵⁴ Reprinted with permission from Hajibabaei et al., Phys. Rev. B 103, 214102 (2021). Copyright 2021 American Physical Society.

Experimentally, Li₇P₃S₁₁ is often prepared in glass-ceramic phase by controlled crystallization and heat treatment of the glass phase. Depending on the heat treatment conditions, ionic conductivities of $3.2 \times 10^{- 3}$ and $1.7 \times 10^{- 2} S / cm$ are reported.^140,141 It was shown that the experimental conductivity is correlated with the degree of crystallization, depending on the heat treatment method.¹⁵⁴ Experimental measurements of the ionic conductivity with techniques such as impedance spectroscopy probe the long-range dynamics of Li, which is influenced by the properties such as grain boundaries and degree of amorphism. Wohlmuth et al.¹⁵⁵ probed the short-range diffusivity of Li in glass-ceramic Li₇P₃S₁₁ using NMR a range of nuclear magnetic resonance (NMR) methods and showed that the activation energy of bulk Li₇P₃S₁₁ is 0.20 eV, which is in agreement with our simulation. It was shown in Sec. IV F that grain boundaries modify the activation energy by $E_{a}^{'} = E_{a} + (E_{g} - E_{a}) / l$ ⁠, where E_g is the energy barrier at the grain boundary and l is the grain size that can often be measured experimentally.

7. Li–X–Y ternary systems

We studied Li diffusivity in Li₇P₃S₁₁ with extensive MD simulations using a ML potential created with the SGPR algorithm. These simulations described the diffusion properties and the melting temperature in agreement with the experimental measurements. We simulated systems with $\sim 700$ atoms for several nanoseconds. A new crystal structure is found for Li₇P₃S₁₁, which is referred to the β phase. The β phase shows much lower ionic conductivity than the α phase. These phases are almost isoenergetic according to exact DFT calculations. Therefore, it is important to identify the β phase experimentally and fine-tune the preparation method in order to avoid this phase. Additionally, we studied diffusivity in the glass phase and showed that the ML potential predicts the crystallization temperature consistently with experimental observations.

In this section, we turn our focus on known Li–X–Y ternary systems, which are available online (materials-project¹¹¹). We considered any crystal that has a considerable concentration of Li to be considered as a solid electrolyte. For any material that qualified certain criteria, a ML potential is generated using the SGPR algorithm and its Li diffusivity properties are calculated. Then, we discuss how the ML models for similar materials can be combined for generating global/universal potentials.

a. Searching algorithm

For all of the Li–X–Y ternary systems, with the exclusion of transition metals, a sparse GP potential was trained using an on-the-fly adaptive sampling algorithm for known crystal structures (Materials Project database) in the materials-project¹¹¹ repository. At the first filtering, the systems with lithium less than 10% of their atoms are dropped. As the next filter, the systems whose electronic bandgap is less than 1 eV are eliminated. After the first two filters, about 300 structures remained for further investigations. Then, on-the-fly ML MD simulations for the potential energy are triggered for these materials. The size of supercells for the ML MD simulations was chosen such that k-spacing $Δ k$ in the reciprocal space satisfied $Δ k < \frac{2 π}{10} Å^{- 1}$ ⁠. Therefore, only Γ-point was chosen as the k-point grid for DFT calculations. The number of atoms in the supercell was in the range of 100–300. If the Li diffusivity of a certain material was estimated to be too small (e.g., less than ∼0.1 Å²/ps at 1000 K), the simulation was stopped for that material. The other condition for continuation of the simulation was stability. If a crystal melted at temperatures higher than ∼600 K, the simulation was terminated unless a very high diffusivity is observed at lower temperatures. The ML MD simulations were carried out in isothermal-isobaric (NPT) ensemble for ∼100 ps. Temperature dependent diffusivity and room temperature ionic conductivity were calculated for this set of materials and listed in Table S1 of Ref. 63.

At the end, 22 promising crystals were detected (Table II). As an exemplary case, 4 materials out of 22 had the same set of elements (Li–P–S) for which the training set consisted of crystal structures known for their relevance to Li-ion conductivity focusing primarily on three materials: γ-Li₃PS₄ with $(1 \times 1 \times 2)$ ⁠, β-Li₃PS₄ with $(2 \times 1 \times 1)$ ⁠, and Li₇P₃S₁₁ with $(2 \times 1 \times 1)$ cells.^138–140 This sampling procedure was discussed in the earlier Subsection V C. As such, the combinatorial approach above opens the door for strategies like “high level parallelism” or “a library of SGPR potentials,” which can be combined by simple cross covariance calculations to describe increasingly complex materials. Such a combinatorial approach is not trivial by other methods.

b. Universal ML interatomic potentials for surveying ternary solid electrolytes

With on-the-fly ML of interatomic potentials using the SGPR algorithm, MD simulations are applied for a survey of Li diffusivity in hundreds of ternary crystals as potential electrolytes for all-solid-state batteries. The models generated for these crystals can be easily combined for creating more general and transferable models, which can potentially be used for simulating new materials without further training. As an example, universal potentials are created for Li–P–S and Li–Ge–S systems by combining the expert models of the crystals that contained the same set of elements. Combinatorial models of different ternary crystals can be directly applied for modeling composite quaternary ones (e.g., Li–Ge–P–S). This hierarchical approach paves the way for modeling large-scale complexity by a combinatorial approach.⁶³

NNs are suitable for large datasets, while kernel-based methods excel in training with smaller data and on-the-fly learning.^{43,44,54,91,101} These methods have found a wide range of application in physical chemistry. Relating to solid electrolytes (SEs), Li₃N is studied using electrostatic spectral neighbor analysis potential (eSNAP) method,¹⁵⁶ Li₃PO₄ using NNs,¹⁵⁷ and Li₇P₃S₁₁ using SGPR.⁵⁴

An important issue for the future of ML IAPs is transferability and systematic model reusability. Currently, models are generated for specific tasks and often discarded afterward. The ability to build on prior work is highly sought-after and imperative for efficient exploration of functional materials. Although gradual improving of the trained models is possible with all ML methods, kernel-based methods such as SGPR offer natural frameworks for libraries of ML models, which can be combined with little effort to describe increasingly complex materials. NNs memorize the data using a set of weights with fixed size and thus have to be re-optimized entirely upon adding new data. In contrast, in kernel-based methods, models are represented by a set of LCEs that can be easily categorized and encapsulated in libraries. These ideas can be explored in the context of Li conducting crystals, which are one of the essential components of all-solid-state storage devices. In an unsupervised fashion, hundreds of ternary crystals are simulated by on-the-fly ML MD and fastest Li conductors are distinguished. Later, we discuss combining the SGPR models for more transferable potentials and a combinatorial approach for modeling new crystals.

Li₇P₃S₁₁ superionic conductor is extensively studied using SGPR. Similar simulations are applied to hundreds of potential SEs. We focus on ternary systems (composed of three elements) in the materialsproject¹¹¹ repository that contains Li. Crystals containing the transition metals are not considered here. Li concentration and bandgap are used as initial filters since a SE should have considerable Li atoms and it should be an electron insulator. About 300 crystals passed the initial filters. A quaternary system of Li₁₀Ge(PS₆)₂¹⁴² is included in the list as an extension of the ternary system. Then NPT-MD simulations with on-the-fly SGPR are performed for a duration of ∼100 ps for each crystal at standard pressure. Three or four temperatures are chosen for MD simulation depending on the stability against melting and the Li diffusion rate.

Figure 16 and Table II show the Li diffusivity in 22 fastest Li conducting crystals. In NPT-MD simulations, the volume expands with increasing temperature. This temperature dependence is ignored in reporting the ionic conductivities in Table II since it only modifies the values by $O (1 %)$ (see Ref. 54). Nevertheless, the effect of NPT ensemble is directly reflected in diffusivity and significantly faster ionic conduction is observed in comparison with NVT ensemble.

FIG. 16.

View large Download slide

Extrapolated ionic conductivity at 300 K (blue) and activation energy (red) for 22 fastest Li conducting crystals. Figure taken from Ref. 63. Reprinted with permission from Hajibabaei and Kim, J. Phys. Chem. Lett. 12, 8115 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

Fast Li diffusion is one of the most important factors for a promising solid electrolyte. Stability against Li metal and wide electrochemical window are the other important factors. Additionally, experimental reports often measure long-range diffusivity, which is affected by grain boundaries, while the values in Table II reflect only the short-range diffusivity. A special technique is required for probing the short-range diffusivity in experimental samples.¹⁵⁵ Among 22 Li conducting SEs reported in Table II, most are sulfides (14), few halides (5), two oxides, and one selenide. Here we briefly discuss the experimental results for selective SEs. The top three fastest ionic conductors belong to the Li₂S-P₂S₅ system. The most stable compound in this system is (γ-/β-) Li₃PS₄. Li₇PS₆ is a lithium argyrodite with orthorhombic/cubic (low/high temperatures) phase where the cubic phase has a higher conductivity. The highest experimental conductivity reported in this system is for Li₇P₃S₁₁, which has a triclinic cell made of PS₄ tetrahedra and P₂S₇ double-tetrahedra.¹⁴⁰ Experimental room temperature ionic conductivities (mS/cm) of 0.16 for Li₃PS_4,¹⁵⁸ 0.11 for Li₇PS_6,¹⁵⁹ and 17 for Li₇P₃S₁₁¹⁴¹ are reported. The latter is one of the highest reported to date for a SE.¹³⁸ As a halide, a conductivity of ∼4 mS/cm (at 330 K) is reported for Li₃InBr₆.¹⁶⁰ Li₄SiO₄ has also been considered as a parent for SE solutions.¹⁵¹ The experimental conductivities are often improved, even by several orders of magnitude,¹⁵⁸ depending on the preparation method. Our simulations suggest that experimental preparation methods could be improved, specially for Li₇PS₆.

8. Universal potentials: Application to quaternary systems

In previous simulations, an “expert” model is generated for the potential energy of each material. At the end, 22 materials are selected, which show significant Li diffusivity. For instance, 4 materials out of 22 have the same set of elements (Li-P-S). It is intriguing to somehow combine these four models to create a universal potential for Li–P–S system. The combinatorial approach based on MOE model opens the door for strategies like “high level parallelism” or “a library of SGPR potentials,” which can be combined by simple cross covariance calculations to describe increasingly complex materials.

As a more simple case, the SGPR models for Li₃PS₄ and Li₄GeS₄ are simply combined and tested on Li₁₀Ge(PS₆)₂ (both configurations in Table II). Due to the presence of different species P and Ge, the combinatorial model cannot be significantly compressed. The combined model of the two ternary crystals described the forces of the quaternary crystal with $R^{2} = 0.961$ ⁠, which is very close to R² of the expert model. This is understandable in view of the similarity of local chemical environments for the cations P and Ge in these crystals (see Fig. 6). In many promising solid electrolytes, they are present as ${(X S_{4})}^{- q}$ poly-anionic moieties with X = P, Ge, etc., in the sea of Li⁺ ions and do not directly interact with each other. By a combination of the ML models for the simpler crystals, we can simulate the composite ones with a minimal training or even without any additional training. Thus, by combining two ternary systems of Li₃PS₄ and Li₄GeS₄, a quaternary system of Li₁₀Ge(PS₆)₂ was studied because the models generated with kernel-based methods such as SGPR can be easily combined with more generality and transferability. Using the previously trained models in the database, it is possible to model new crystals without further training. This approach is indeed demonstrated by creating universal potentials for two ternary systems. Combinations of two ternary crystals resulted in a realistically accurate potential for a quaternary crystal. With the introduced algorithm, the models generated in this study can be reused, with limited or no training at all, to describe more complex system (such as quaternary electrolytes) by a combinatorial approach.

9. Summary for SGPR simulations of solid electrolytes

The SGPR-based ML formalism is applied to the studies of ionic diffusion in dozens of solid electrolytes. Specifically for Li₇P₃S₁₁, these simulations described the diffusion properties and the melting temperature in agreement with experiments. We simulated systems with $\sim 700$ atoms for several nanoseconds, which is unimaginable with direct ab initio molecular dynamics (AIMD). Although we could easily simulate larger systems because the computational cost for calculating the energy and forces is proportional to the number of atoms N, unlike DFT that scales as N³ or $N^{2} log N$ ⁠. For $\sim 700$ atoms, using 100 CPUs, it took less than 1 s to predict the ab initio potential energy and forces using the SGPR potential which was faster than direct DFT by a factor of $O (10^{4})$ ⁠. A new crystal structure is found for Li₇P₃S₁₁, which is referred to as the β phase. The β phase shows much lower ionic conductivity than α phase. These phases are almost isoenergetic according to exact DFT calculations. Therefore, it is important to identify the β phase experimentally and fine-tune the preparation method in order to avoid this phase. Additionally, we studied diffusivity in the glass phase and showed that the ML potential predicts the crystallization temperature consistently with the experimental observations.

We showed that SGPR models for materials with the same set of elements can be easily merged into a single universal potential. This allows high level parallelism in generating ML potentials, which is not easily possible with other ML algorithms. That is, instead of training a potential sequentially in different domains, we can train potentials for all of the domains in parallel and merge them into a single global potential at the end. This is a critical step forward for applications of the ML potentials at increasingly wider domains and truly complex materials.

B. Lithium ion batteries

1. Al-doping in Li-ion-battery cathode materials

Al-doping has been employed in the cathode materials of lithium-ion batteries (LIBs) to alter the reaction of negatively charged ions in order to prevent declining in energy storing ability or dropping in electrical potential. The chemistry involving the transfer of negative charge in cathodes with high lithium content has opened avenues for creating cathode materials with greater energy capacity. Yet, the prolonged reaction of these negatively charged ions leads to a reduction in the battery's storage capacity and a decrease in voltage as the battery is used repeatedly. This occurs due to irreversible changes in the oxygen molecules, destabilizing the structure. To enhance the durability and steadiness of layered oxides and lithium-rich variants such as lithium nickel–cobalt–aluminum oxide (LNCAO or abbreviated as NCA) and lithium nickel–cobalt–manganese oxide (LNCMO or abbreviated as NCM), Al-doping has been intermittently applied. However, prior studies on Al-doping have mostly focused on NCA or NCM, which predominantly contain TM elements from the 3d group. The chemical behavior differs significantly in the case where the elements from the 4d or TM groups are utilized, leading to distinct characteristics. Cathodes comprising solely 3d TMs exhibit unique structural arrangements. In contrast, addition of more 4d or 5d transition metals into oxide-based cathodes introduces novel chemical properties. Due to the intricate behavior of these materials, SGPR-based DFT potentials/forces were employed to identify their most stable structures, presenting a new approach to designing fully reversible layered cathode materials. By employing this method, the impact of Ni and Al-doping was investigated on the decline in capacity and voltage of the high-capacity LIB cathode materials containing 4d-elements. An exhaustive and dependable analysis of this complex system was conducted using SGPR.⁶⁵ The instability of LRO (Li₂RuO₃; Li_1.33Ru_0.67O₂) was significantly enhanced by Ni-doping (LRNO: Li_1.22Ru_0.61Ni_0.17O₂). Further, stability was achieved in this LRNO through additional Al-doping. These Al-doped LRNO (Al-LRNO: Li_1.22Ru_0.61Ni_0.11Al_0.06O₂) materials exhibit exceptional stability and high capacity. The structures of LRNO and Al-LRNO were investigated with SGPR-based DFT potentials by incorporating a diverse genetic algorithm based on the known LNO structure. The effects of N-doping and Al-doping on the materials' stability under high voltage conditions were investigated.

2. DFT calculations

To generate training test sets for SGPR potential energies, spin-polarized DFT calculations were conducted to derive potential energies for both discharged and charged phases of LRNO and Al-LRNO configurations by using a PAW plane wave with an energy cutoff of 600 eV. The method involved PBE+U functional with Tkatcheko–Schaffler (TS)¹⁶¹ dispersion correction, where U = 7.5 eV, J = 0.6 eV for Ru and U = 7.7 eV, J = 1.0 eV for Ni. The structures were fully relaxed until the interatomic forces were smaller than 0.03 eV/Å.

3. Structure prediction

To identify the best fitting arrangements of cations in the LRNO and Al-LRNO systems, the structures with the lowest energy were searched by altering atomic positions within the unit-cell formula to mirror the effect of doping. A modified genetic algorithm was employed starting with a handful of randomly generated parent structures. These parents represented changes in atomic positions, swapping three Li/Ru-sites within the 3 × 2 × 1 Li₂RuO₃ supercell with either three Ni atoms or a combination of two Ni and one Al atom, thereby forming LRNO or Al-LRNO reference systems, respectively. At each step, numerous child structures were randomly generated from each parent, and the lowest energy structures were selected from the combined parent and child generation for the next iteration. The search process continued until minimal changes occurred within the parent structures over a few steps, while conducting multiple independent searches. To calculate the energy of each structure, the ML potential was generated on-the-fly-by using the SGPR methodology, employing a variation of SOAPs—a many-body descriptor—as the core kernel. During each energy calculation step, if the ML model fell short, precise DFT calculations were performed, updating the model. An accuracy of 3 meV per atom was achieved by retaining the structures that could potentially possess the minimum energy within this margin for exact DFT calculations. Geometry optimization was performed to determine the lowest energy structure. Upon identifying stable cation ordering in Al-LRNO, a Li vacancy was introduced to balance additional charge, and its placement was determined by defect formation energy. To forecast the charging phase, each state of charge was handled separately, thus reducing the number of sampled geometries. Commencing with pre-established cation ordering, ∼500 random delithiated geometries were generated for each model, with their structure energies predicted by ML potential. The SGPR model recognized unique environments from around 11 samples among prepared structures, learning atomic potential from DFT results. Again, the model identified the Li ordering with the lowest energy, validated subsequently by DFT. The trained model underwent validation by comparing predicted and calculated formation energies of six randomly selected non-learning samples. For the training/test sets, each model exhibited R² of 0.96/0.92, RMSE of 10.2/10.4 meV, and MAE of 8.2/6.4 meV (Fig. 17).

FIG. 17.

View large Download slide

Predicted formation energy (⁠ $E_{form}^{M L}$ ⁠) vs DFT formation energy (⁠ $E_{form}^{DFT}$ ⁠) of train set and test set at each delithiation. (a) 50% delithiated LRNO, (b) 75% delithiated LRNO, (c) 50% delithiated Al-LRNO, and (d) 75% delithiated Al-LRNO.⁶⁵ Figure taken from Ref. 65. Reprinted with permission from Ha et al., Adv. Energy Mater. 12, 2201497 (2022). Copyright 2022 Authors, licensed under a Creative Commons Attribution (CC BY) license.

4. Structures of LRNO and Al-LRNO

The structure features a honeycomb-like Li-transition metal layer reminiscent of LRO, showcasing an edge-shared section composed of RuO₆-NiO₆ or RuO₆-AlO₆ configurations [Fig. 18(a)]. A resemblance is observed in their respective calculated x-ray diffraction (XRD) spectra [Fig. 18(b)]. When comparing the optimized LRO, predicted LRNO, and experimental XRD results, no significant visual disparities were noted.^162,163 This alignment significantly bolsters the reliability of the SGPR-based ML-driven simulations rooted in first principles for the structures of LRO, LRNO, and Al-LRNO. A notable observation is that in LRNO, Ni tends to form clusters, whereas in Al-LRNO, Ni and Al coalesce into alloy-type clusters, predominantly bonding with each other rather than exclusively with Ru. This inclination is supported by the cluster expansion coefficients (energies) derived from pair atoms, obtained through the least absolute shrinkage and selection operator (LASSO) regression model.¹⁶⁴ This intriguing pattern can be comprehended through the concept of immiscibility or miscibility between Al, Ni, and Ru.^166,167 DFT calculations on the alloying or mixing energies of these metals suggest that Ni and Ru tend to be immiscible, leading to the clustering of Ni atoms adjacent to Ru in LRNO. Meanwhile, the calculations indicate that Al–Ni mixing is marginally more favorable than Al–Ru mixing, resulting in a preference for Al atoms to bond primarily with Ni atoms, especially when a small quantity of Al atoms is present.

FIG. 18.

View large Download slide

(a) Simulated atomic structures (top/side views) of LRO, LRNO, and Al-LRNO crystals generated by SGPR-based ML methods along with genetic algorithm approach and (b) their XRD patterns. The XRD patterns are calculated using λ = 1.541 Å. The arrangements of various atoms are shown in color code: green Li; skyblue Ru; orange Ni; deep-blue Al; red O.⁶⁵ Figure taken from Ref. 65. Reprinted with permission from Ha et al., Adv. Energy Mater. 12, 2201497 (2022). Copyright 2022 Authors, licensed under a Creative Commons Attribution (CC BY) license.

5. Stability fading of Al-LRNO due to delithiation at high-voltages

The delithiated structures of LRNO and Al-LRNO are investigated to scrutinize the impact of Al-doping on stability, particularly in high-voltage environments. Initially, the fully (100%) delithiated structures of LRNO and Al-LRNO are optimized. Delithiation in LRO triggers significant distortion, indicated by a short O–O distance of approximately 2.39 Å, resulting from anionic activity.¹⁶⁷ This leads to an essential de-coordination of oxygen atoms, inducing a structural distortion to alleviate strain energy through a reductive coupling mechanism. In LRO, oxygen prompts adjacent oxygen atoms to approach Ru along the O–O axis, causing O–O dimerization and eventual O₂ release, contributing to poor anionic cycling stability, as highlighted in prior studies. However, in LRNO after Ni substitution, the distortion is somewhat alleviated, evident from a larger O–O distance (∼2.42–2.58 Å) compared to LRO. This suggests suppression of oxygen oxidation. Subsequent Al-doping in LRNO relaxes the delithiated structure further, resulting in a significantly increased O–O distance (∼2.50–2.62 Å). The O atoms attached to Al, exhibiting strong ionic covalent characteristics, hinder O–O dimerization, thereby suppressing O₂ release. Additionally, the O defect formation energy was calculated to comprehend the suppression of structural distortion due to O₂ release. The calculations reveal a higher O defect formation energy in Al-LRNO > LRNO > LRO at any delithiated state. By recognizing that 100% delithiation during charging is unrealistic, the 50% and 75% delithiated structures of LRNO and Al-LRNO were examined to determine their formation energies (meV/f.u., where f.u. is formula-unit) using SGPR model-based DFT-level potential. The results confirm the higher stability of delithiated Al-LRNO structures over the corresponding delithiated LRNO structures, evidenced by more negative formation energies [Fig. 19(a)]. Despite similar destabilization from structural reconstruction in both Al-LRNO and LRNO, the heightened stability in delithiated Al-LRNO stems from Al's stronger reduction effect compared to Ni, as validated by the Bader charges. Specifically, the Bader charges (in au) of O around Ni in LRNO for 50% and 75% delithiated phases are −1.01 and −0.87, respectively. Conversely, around Al in Al-LRNO for 50% and 75% phases, the values become more negative, registering at −1.23 and −1.12, respectively. These findings imply Al effectively reinforces oxygen binding energy in Al-LRNO, leading to greater suppression of oxygen evolution. Furthermore, a notable difference is observed in delithiation around the Al atom compared to LRNO, where delithiation is significantly curtailed [Figs. 19(b)–19(g)]. At highly delithiated states, Li-ions farther from Al tend to diffuse out first, whereas no specific Li position preference is noted around Ni.

FIG. 19.

View large Download slide

Machine learning predicted structures and formation energy of LRO, LRNO, and Al-LRNO at discharged and charged states. (a) Formation energy (meV/f.u.) with delithiation process in LRNO (red line) and Al-LRNO (blue line). The most stable geometry (filled circle) is found among sampled configurations (empty triangle) by the sparse Gaussian process regression (SGPR) model-based on-the-fly potential learning algorithm. (b) and (c) The SGPR model-based machine learning predicted pristine structures of (b) LRNO and (c) Al-LRNO at the discharged phase. (d)–(g) Charging phases (50% and 75%) (Li: green, Ru: skyblue, Ni: orange, Al: blue, O: red).⁶⁵ Figure taken from Ref. 65. Reprinted with permission from Ha et al., Adv. Energy Mater. 12, 2201497 (2022). Copyright 2022 Authors, licensed under a Creative Commons Attribution (CC BY) license.

6. Enhanced diffusion by Al-doping

Delving into the diffusion coefficient via AIMD study poses challenges due to the fluctuating lithium concentration during the charging and discharging cycles. Particularly, at room temperature, the diffusion of Li atoms belongs to extremely rare events. To expedite computations, the diffusion was simulated at elevated temperatures of 700, 800, 900, 1000, and 1200 K. To find the impact of Al-doping, LRNO and Al-LRNO structures were examined by employing (2 × 1 × 1) cells with a 50% Li concentration, comprising 122 and 120 atoms in the simulation cells, respectively. The SGPR-based ML DFT potential was trained dynamically during the simulations. MD simulations in the NVT ensemble with a Nosé–Hoover thermostat were conducted for 0.5 ns (except for 1 ns at 700 K), using a time step of 2 fs. Each atom's mean squared displacement (MSD) from its initial position surpassed 50 Å within the simulation timeframe. Typically, the temperature dependency of the diffusion constant D is approximated using the Arrhenius relation: $D = D_{o} exp (- E_{a} / k_{b} T)$ ⁠, where E_a represents the activation energy. The calculated activation energies for LRNO and Al-LRNO with 50% Li concentration stand at 0.48 and 0.40 eV, respectively (Fig. 20). Simulations based on the SGPR ML DFT potential project a higher diffusion coefficient for the Al-doped material compared to the pristine one, suggesting an approximate doubling in Li⁺ mobility within the cathode material, especially around ∼700 °C, where Li⁺ diffusion favors Al-LRNO over LRNO.

FIG. 20.

View large Download slide

Diffusion coefficient of Li ions (with 50% concentration) in LRNO (blue bullets) and Al-LRNO (red bullets) calculated with SGPR-based machine-learned DFT potential MD simulations.⁶⁵ Figure taken from Ref. 65. Reprinted with permission from Ha et al., Adv. Energy Mater. 12, 2201497 (2022). Copyright 2022 Authors, licensed under a Creative Commons Attribution (CC BY) license.

7. Electronic origin of stabilization and Discussion

In the realm of cathode materials containing 4d-elements, the introduction of Al through doping has shown remarkable efficacy in mitigating capacity and voltage deterioration. This intervention not only addresses persistent challenges related to capacity fading but also curtails the decline in average discharge voltage, thereby paving the way for high-performance LIBs. Al-doping fosters the formation of Ni–Al bonding, bolstering the characteristics of bonding orbitals within Al–O bonds. This fortified Al–O bonding acts as a deterrent to oxygen oxidation, significantly enhancing structural stability and alleviating safety concerns. The capability of Al-doping to suppress capacity fading and voltage decay is pivotal in the development of stable, reversible layered cathode materials. Furthermore, the reported material possesses the potential to serve as a protective coating over cathode materials, substantially enhancing their electrochemical performance by imbuing them with heightened stability. The exceptional performance of this material can be attributed to several factors. The bonding orbital characteristics of Al–O, situated deep below the Fermi level in Al-doped material, effectively restrain excessive oxygen oxidation even during profound delithiation. This phenomenon stabilizes crystal structures by notably suppressing oxygen defect formation compared to LRNO. Al-LRNO demonstrates higher structural stability within the high voltage region due to the less oxidized state of Ni, a characteristic observed in this Al-doped material. In essence, the strategic utilization of Al-doping stands as a promising avenue for addressing key challenges in LIBs, fostering improved stability and performance.

Notably, this SGPR ML-based computational approach underscores the effective suppression of cationic migration within these high-voltage redox systems. Al-doping emerged as a highly impactful strategy, influencing the bonding dynamics to stabilize TM octahedrons within the high-voltage redox framework. The Al-doping not only enhances diffusion rates but also reduces Li-ion charge transfer resistance, bolstering the overall performance. Drawing from these fundamental insights and preceding experimental results, conventional cathode materials such as NCMA could potentially be harnessed for long-term development in high-energy-density batteries. The analogical relationship between Mn and Ru in the periodic table holds promise for this purpose. The exceptional attributes of Al-LRNO, including its high capacity, diminished voltage hysteresis, and remarkable structural integrity, owe credit to the induced Al–O bonding orbital characteristics. The ionic covalent nature of Al–O bonding orbitals, situated significantly below the Fermi level within the doped material, plays a pivotal role in restraining excessive oxygen oxidation, even under deep delithiation conditions, effectively stabilizing the crystal structure. We posit that the reported material holds potential as a protective coating for commercial cathode materials, fostering high-temperature stability and substantially augmenting their electrochemical performance. This strategic use of the material may pave the way for advancements in battery technology.

8. Sodium-substituted disordered rock salt cathode: Reduced energy barrier for oxygen evolution

Cation-disordered rock salt (DRX) cathodes stand as promising candidates for future high-energy-density materials, surpassing conventional layered cathodes in LIB technology. An intriguing discovery emerged when Na-doped DRX exhibited efficient electrocatalytic properties toward the OER. This OER electrocatalyst demonstrated a current density of 10 mA cm⁻² at an overpotential (η) of 270 mV, displaying a Tafel slope of 67.5 mV dec⁻¹ and impressive long-term stability surpassing 5.5 days—surpassing the benchmark IrO₂ (with η = 330 mV and Tafel slope = 74.8 mV dec⁻¹). This intriguing electrochemical behavior finds strong support both experimentally and through ML-based searches for the minimum energy structure. Na-doping notably facilitates O₂ evolution by reducing the oxygen binding energy (OBE) on the surface, showing a decrease from 6.51 eV in pristine DRX to 5.45 eV with Na-doping. To comprehend Na's role, the scenario where 1/16 of Li atoms in Li₂RuO₃ were substituted by Na was examined. Simulation-based relaxations considering positions and cell structures demonstrated that the system's potential energy is minimized when the Na atom occupies the Li layer rather than the Ru layer, due to the compatibility between Na and Li. Subsequent discharging/delithiation was simulated using AIMD. SGPR MD was performed to construct an on-the-fly potential energy surface of the spin-polarized PBE+U calculations.⁶⁶ The ML model accurately predicted the DFT forces with a RMSE of ∼0.2 eV, which can be considered favorable at high temperatures. Simulations were conducted with a time step of 2 fs in the isothermal isobaric ensemble (NPT) to monitor volume and shape changes during discharging at an elevated temperature of 700 K for enhanced dynamics. Two supercells were created—one with 32 undoped Li atoms and another with two Li atoms replaced by Na atoms in their minimum energy positions. Both were simulated using MLMD to discern differences. During these simulations, the delithiation potential was calculated based on energy changes and volume variations in the crystal structure. The findings highlighted that while Na's presence did not notably affect the potential, it played a pivotal role in controlling crystal volume. At lower delithiation stages (<50%), the Na-doped structure exhibited slightly higher volume than the undoped crystal. However, in higher delithiation percentages (50%–80%), the presence of Na prevented excessive expansion along the stacking direction of the Ru layers, thereby enhancing the system's mechanical stability (Fig. 21).

FIG. 21.

Effect of Na doping in DRX. (Left) Supercell of Na-doped Li2RuO3 used for delithiation (discharging) simulations with SGPR based PBE+D3 first-principles MD simulation. After every 5 ps a Li atom was removed from system. (Right) Spacing of Ru layers and delithiation potential averaged over every 5 ps. All calculations were spin-polarized and the maximum k-spacing of |Δk|<π/10Å−1 was applied for generating the k-point grid. A kinetic energy cutoff of 500 eV is applied. Due to the large size of unit cells, the Γ-centered k-Point grid of (1, 2, 3) was chosen. The (U; J) values of (4:5; 0:6) for Ru and (4:7; 1:0) for Ni were used for PBE+U calculations.66

View large Download slide

Effect of Na doping in DRX. (Left) Supercell of Na-doped Li₂RuO₃ used for delithiation (discharging) simulations with SGPR based PBE+D3 first-principles MD simulation. After every 5 ps a Li atom was removed from system. (Right) Spacing of Ru layers and delithiation potential averaged over every 5 ps. All calculations were spin-polarized and the maximum k-spacing of $| Δ k | < π / 10 Å^{- 1}$ was applied for generating the k-point grid. A kinetic energy cutoff of 500 eV is applied. Due to the large size of unit cells, the Γ-centered k-Point grid of (1, 2, 3) was chosen. The (U; J) values of (4:5; 0:6) for Ru and (4:7; 1:0) for Ni were used for PBE+U calculations.⁶⁶

C. Electrocatalysts

1. High-performing atomic electrocatalysts

There is an urgent need for developing efficient and environmentally friendly electrocatalysts for various electrochemical reactions including hydrogen/oxygen evolution/reduction, carbon dioxide reduction, ammonia reduction, hydroxide evolution, and chlorine evolution. Traditionally, these catalysts relied on expensive metals like platinum, ruthenium, or iridium. Recent efforts focus on using minimal amounts of these costly metals or finding cheaper alternatives, embedding them into graphene or other 2D materials to maximize surface area and catalytic activity. However, while this embedding enhances the active surface area, it can decrease electrical conductivity, impacting the catalyst's efficiency. To optimize these catalysts, researchers turn to high-throughput computer screening, analyzing different metal configurations onto/into 2D materials and their impact on stability, electronic properties, and electrochemical reduction/evolution activities. Their findings highlight specific configurations as displaying superior durability and high conductivity for efficient water splitting, crucial for clean energy applications like hydrogen-powered vehicles and electricity generation. The study also explores ML approaches and AIMD simulations to characterize high-performance catalysts.^168–170 Detailed examinations have unveiled that single atom catalysts (SACs)^161,172 follow distinct catalytic pathways compared to their nanocluster counterparts, resulting in differing selectivity in catalytic reactions. This means that by precisely adjusting the coordination environment of active sites, specific catalytic reaction pathways and product preferences can be achieved. Notably, many carbon-based materials have proven to be effective hosts for anchoring single transition metal atoms, serving as a favorable environment for these catalysts.⁶⁷

Lately, there has been a surge in interest toward atomically dispersed catalysts that feature isolated transition metal atoms supported on a variety of materials such as metal oxides, borides, carbides, nitrides, graphene, graphitic carbon nitride (g-C₃N₄), hexagonal-boron nitride (h-BN), black phosphorous (BP), layered double hydroxides (LDHs), and transition metal dichalcogenides (TMDs), MXenes, MBenes.¹⁷³ These catalysts have garnered attention due to their notably high utilization of atoms, exposed active sites, customizable electronic structure properties, and robust interaction between the metal and the supporting material. Known as single atom catalysts (SACs),^161,172 embedding atoms catalysts onto or into 2D materials are being widely explored and recognized as promising candidates for a wide-ranging catalytic reactions such as hydrogen evolution (HER),^{165,174–177} oxygen evolution/reduction (OER/ORR),^{166,178–180} nitrogen reduction (NRR),^{173,181–184} CO/CO₂ reduction (CORR/CO₂RR),⁶⁹ and chlorine evolution reaction (ChER).¹⁸⁵ Ultimately, this research contributes to a clearer understanding of these catalysts, particularly for critical applications in clean energy environments like hydrogen fuel cells and electric power generation.

2. HER/OER/ORR

Tiny yet mighty, nanoparticle (NP)-sized noble electrocatalysts have been harnessed for a range of electrochemical reactions, particularly in advancing an environmentally friendly hydrogen-based economy like water splitting. Researchers have leveraged minuscule quantities of single atoms (SAs) to amplify the active surface area, fine-tuning catalytic abilities by situating these atoms within defect sites on nitrogen-doped graphene (G_N). Developing highly effective catalysts that serve multiple purpose like HER, OER, and ORR has become a focal point in advancing clean and sustainable energy technologies.¹⁷¹ Focusing on these reactions, top-tier performances from 3d-5d TM SACs have been unveiled by employing DFT and ML. The investigation delved into the stability and activity of TM–G_N, scrutinizing structure/coordination, formation energy, structural/electrochemical stability, electronic traits, conductivity, and reaction mechanisms. Among the array of –N_nC_m configurations, the –N₂C₂ compositions stand out for their ease of formation, showcasing superior electrochemical catalytic prowess and prolonged durability, in particular, for HER, notably avoiding aggregation or dissolution, compared to the extensively studied –C₄/C₃ and –N₄/N₃ configurations.⁶⁵

Introducing non-noble TMs such as Fe/Ni/Cu onto the surface of MXenes as individual SAs has demonstrated altered mechanical and electrical traits, albeit rarely observed. Through the use of DFT simulations, various combinations of MXenes are assessed for their stability with intermediates and adsorption free energies to gauge their efficacy as catalysts for these reactions. The majority of the considered MXenes remain thermodynamically stable upon doping with Ni/Cu/Fe, showcasing optimal performance for OER/ORR/HER, respectively. The ML potential was generated using on-the-fly sampling and SGPR algorithm.⁷⁰ The kinetic stability of the material was confirmed by using SGPR ML potential for the catalyst. Surprisingly minimal overpotentials observed in these reactions highlight these catalysts' appeal for water splitting (overall potential of 0.48 V) and metal–air batteries (overall potential of 0.70 V), surpassing even commercially utilized IrO₂(111) and Pt(111) for the OER and HER/ORR, respectively.^70,186 The strong correlations observed between intermediates (scaling relations) and their interactions with the catalytic surface predominantly dictate their performance. Doped metals not only serve as active sites but also play a crucial role as conduits for charge transfer between the surface and adsorbate. The notably low reaction barriers of Ni–Cr₂ScC₂O₂ for all three reactions, coupled with robust reaction kinetics, position it as an efficient catalyst serving three functions simultaneously.

Here, we show heterogeneous catalysis dynamics of H/Pt reported by Vandermause et al.⁴⁵ They validated a ML sparse Gaussian process (SGP) force field model for simulating Pt, H₂, and H/Pt systems, comparing its performance with DFT and ReaxFF^187,188 calculations. The SGP model was trained on energy, forces, and stress data from MD simulations and tested on various properties not included in the training set. The SGP model predicted the lattice constant within 0.1% and the bulk modulus and elastic constants within 6% of DFT values. It performed significantly better than ReaxFF, particularly for the C₄₄ elastic constant. The model accurately predicted energies and properties for structures not encountered during training, showing low uncertainties and excellent agreement with DFT, especially for bulk Pt equilibrium volume and H₂ bond lengths. It generalized well for different Pt surface facets. The SGP model showed higher fidelity to DFT, avoiding non-physical behaviors like surface sublimation of Pt hydride species observed in ReaxFF. The model's predictions of H₂ dissociative adsorption and atomic H diffusion on Pt(111) closely matched DFT results. MD simulations using SGP yielded an apparent activation energy and diffusion prefactor consistent with experimental values. Overall, the SGP model demonstrated superior accuracy and reliability in comparison to ReaxFF and showed promising generalization capabilities for unseen configurations (Fig. 22).

FIG. 22.

View large Download slide

Validation of the H/Pt SGP force field. In each plot, shaded regions indicate 99% confidence regions of the SGP model. (a) Energy vs volume of bulk Pt. (b) Energy vs H–H bond distance of a single H₂ molecule. (c) Energy vs intermolecular distance of two H₂ molecules oriented perpendicular (blue) and parallel (green) to each other. (d) Surface energies of Pt for (111) (110) (100), and (211) facets. (e) H adsorption energies at different binding sites and Pt facets (copyright npg Nature).⁴⁵ Figure taken from Ref. 45. Vandermause et al., Nat. Commun. 13, 5183 (2022). Copyright 2022 Authors, licensed under a Creative Commons Attribution (CC BY) license.

3. CO₂RR

Carbon capture, storage, and conversion technologies are seen as viable strategies to tackle excessive CO₂ emissions. The electrocatalytic reduction of CO₂ into carbon-based fuels and chemicals using renewable electricity is a promising method to mitigate CO₂ levels in the atmosphere and address energy concerns. While carbon monoxide (CO) and formic acid (HCOOH) are common products of CO₂ reduction, further transformations can yield valuable outputs like methane, methanol, ethylene, and ethanol. However, CO₂'s chemical inertness leads to high overpotentials and low selectivity for desired products, making the process challenging. High voltages may trigger hydrogen evolution and hinder CO₂ reduction efficiency. While carbon-based materials effectively host TM SAs, stability during catalysis remains an issue. Efforts are focused on identifying non-carbon 2D materials with strong metal-support interaction and structural tunability. For instance, ultrathin gallium nitride (GaN) monolayers have been synthesized as promising hosts for TM SAs, showcasing exceptional electronic and optical properties. GaN-based SACs have demonstrated efficacy in various catalytic reactions such as nitrogen reduction, oxygen evolution/reduction, and more, showcasing potential multifunctionality through precise coordination modulation of active sites. This advancement highlights the potential of GaN SACs in diverse catalytic reactions by precisely adjusting the coordination environment of active sites.⁶⁹

AIMD simulations based on SGPR ML potential at 500 K using the Nose–Hoover thermostat method tested the kinetic stability of proposed structures. The representative examples of (a) Ir-SAC and (b) Os-SAC also showed the dynamic stability of these structures (Fig. 23).

FIG. 23.

View large Download slide

Structural stabilities of proposed TM-SACs: Representative structure of GaN monolayer surface with TM embedded at Ga-SV defect site. Snapshots of the ab initio molecular dynamics simulation at different time. The temperature is controlled at 500 K. The structural stabilities of the representative Ir-SAC and Os-SAC are dynamically stable (Ga: light green, Os: light pink, N: blue).⁶⁹ Figure adapted from Ref. 69. Reproduced with permission from Umer et al., J. Mater. Chem. A 10, 24280–24289 (2022). Copyright 2022 Royal Society of Chemistry.

4. NRR

Nitrogen (N₂) is abundant but inert, requiring transformation into bio-available forms like ammonia (NH₃) and other nitrogen compounds through a process called nitrogen fixation. This process is crucial for sustaining life, impacting agriculture, soil fertility, ecosystem health, and even potential energy sources. NH₃, a key product of nitrogen fixation, is vital in fertilizer production, pharmaceuticals, and energy storage due to its high energy density and hydrogen content. As the global population grows, the demand for NH₃ escalates, necessitating cost-effective production methods that are environmentally friendly. Current industrial methods, like the Haber–Bosch process, though pivotal for food production, contribute significantly to greenhouse gas emissions. Efforts are underway to explore alternatives to the energy-intensive and environmentally impactful industrial processes. Electrochemical nitrogen reduction (NRR) shows promise as a cleaner, more sustainable method for NH₃ synthesis, drawing inspiration from biological nitrogenase enzymes. However, challenges persist in finding efficient catalysts that can enhance the selectivity and yield of NH₃ production while competing with hydrogen evolution reactions. Researchers are investigating diverse materials, including transition metals, borides, sulfides, oxides, single-atom catalysts, MBenes, and defect-engineered substrates, to improve NRR efficiency. Various studies aim to elucidate mechanisms, control parameters, and catalyst performance using experimental, theoretical, and computational approaches and to explore innovative solutions to improve electrochemical NRR performance in the coming years.

SACs for effective nitrogen reduction reaction (NRR) using various transition metals were investigated as implanted on boron-arsenide (BAs), boron-phosphide (BP), and boron-antimony (BSb) using DFT. The kinetic stability of the material was investigated by using simulations based on SGPR machine-learned potential for the catalyst (W-BAs) with the lowest potential determining step (PDS). W-BAs showed high catalytic activity and excellent selectivity with an insignificant barrier of only 0.05 eV along the distal pathway and a surmountable kinetic barrier of 0.34 eV. The W-BSb and Mo-BSb exhibit high performances with limiting potentials of −0.19 and −0.34 V, respectively.

To check the dynamic stability of single and double vacancy systems, SGPR based ML DFT potential/force MD simulations were performed at 500 K. The simulations indicated the dynamical stability of (a) Pd@2DCP, (b) Fe@h-BN, and (c) Fe@h-B₂N₂ structures (Fig. 24),^181–183 which promised the stability of mono- and dual-type non-metal doped TM-SACs which can be utilized for NRR.

FIG. 24.

View large Download slide

Temperature and energy profiles from 6 ps ab initio molecular dynamics simulations under 500 K with the time step of 2 fs. Top and side views of (a) Pd@2DCP, (b) Fe@h-BN, and (c) Fe@h-B₂N₂ structures are presented. Color code: Pd, light-gray; Fe, golden-brown; B, light pink; N, blue; C, gray; O, red; H, cyan).⁶⁸ Figure adapted from Ref. 68. Reproduced with permission from Umer et al., J. Mater. Chem. A 10, 6679–6689 (2022). Copyright 2022 Royal Society of Chemistry.

5. ClER

The chlor-alkali process stands as one of the fundamental chemical industries globally, with an annual chlorine production surpassing 75 × 10⁶ tons. The production of each ton of chlorine demands approximately 2500 kW $\cdot$ h of electricity, contributing to an overall requirement of around 200 TW $\cdot$ h annually for the chlor-alkali sector alone. This consumption constitutes roughly 1% of the world's total energy production. Given these substantial energy demands, it becomes paramount to explore novel, energy-efficient, and eco-friendly methods for synthesizing Cl₂ at a lower cost. Furthermore, achieving selectivity in the oxygen evolution reaction over the chlorine evolution reaction is crucial for industries like pharmaceuticals, cosmetics, and computing. A promising and cost-effective technique to produce pure Cl₂ involves the use of a membrane-separated device within the chlor-alkali process. Presently, the electrocatalysts based on RuO₂ and IrO₂ dominate chlorine evolution reaction (ClER). However, these metal oxides suffer from inadequate catalytic activity, selectivity, and durability due to concurrent oxygen generation during Cl₂ production. For highly efficient ClER, a catalyst capable of producing Cl₂ with minimal overpotential in the chlor-alkali process is imperative. In this regard, for more efficient electrocatalysts, ML potential/force-fields simulations were carried out. The optimal configurations of the complex systems were obtained by simulation/optimization process utilizing the SGPR based ML DFT potential. The free energy profiles of ClERs for Pt–C₄, Pt–N₄, Pt–C₂N₂, Pt–C₂ $N_{2}^{- 2}$ ⁠, and Pt–C₂ $N_{2}^{- 3}$ are shown in Fig. 25. In this context, the Pt–C₂N₂ moiety steered electrocatalysts (referred to as Pt-1) exhibit promising ClER features. Comparative studies demonstrate that Pt-1 displays an onset potential of only 5 mV at 80 °C, in line with industrial Cl₂ production temperatures. At a current density of 10 mA/cm², it exhibits a minimal overpotential (η) of only 39 mV, contrasting with ∼100 mV for the best industrial electrocatalyst. These results align closely with the DFT study conducted. Pt-1 showcases superior selectivity in Cl₂ production (up to 100%, compared to 90% for industrial electrocatalysts), exceptionally high mass catalytic activity (∼2.12 A/mg_Pt against 1.49 × 10⁻⁵ A/mg_Ir−Ru at 70 mV overpotential), and long-term durability in acidic electrolytes. Notably, the reported ClER mass activity for Pt-1 in acidic media surpasses that achieved by any other catalyst to date. The efficient Cl₂ generation capability of Pt–C₂N₂ at remarkably low applied voltages in acidic conditions signifies a groundbreaking advancement toward eco-friendly and economical pure Cl₂ production, potentially saving substantial energy. The as-prepared Pt-N₂C₂ catalyst exhibits superior electrochemical performance for chlorine production in acidic electrolytes. These catalysts can be industrially synthesized at low temperatures, demonstrating exceptional activity, durability, and selectivity (100%)—outperforming commercial and membrane catalysts commonly adopted in industry. The industrial electrodes (IE) comprise substantial amounts of precious metals (⁠ $0.8 - 4$ mg/cm²). Their deactivation arises from weak adhesion of the coating to the substrate, leading to the separation of electrocatalyst species from the electrode surface. Additionally, oxygen penetration through cracks during electrolysis causes electrode passivation. This process reduces the electroactive sites due to the oxidation of RuO₂ and IrO₂ to RuO₄ and IrO₃ at high potentials. However, the Pt-C₂N₂ species in our electrode material circumvent these issues, resulting in heightened stability (Fig. 25).

FIG. 25.

DFT Free energy profiles of ClERs for SGPR based ML-DFT optimized models of Pt–C4, Pt–N4, Pt–C2N2, and PtO2 (110). Pt–C2N2 shows 0.05 eV (and −0.02 eV at 80 °C) at U = Ueq(ClER) in acid condition. (H: white, C: gray, N: deep-blue, O: red, Cl: green, Pt: light-blue).185 Figure adapted from Ref. 191. Reproduced with permission from Ha et al., Small 19, 2300240 (2023). Copyright 2023 Authors, licensed under a Creative Commons Attribution (CC BY) license.

View large Download slide

DFT Free energy profiles of ClERs for SGPR based ML-DFT optimized models of Pt–C₄, Pt–N₄, Pt–C₂N₂, and PtO₂ (110). Pt–C₂N₂ shows 0.05 eV (and −0.02 eV at 80 °C) at U = U_eq(ClER) in acid condition. (H: white, C: gray, N: deep-blue, O: red, Cl: green, Pt: light-blue).¹⁸⁵ Figure adapted from Ref. 191. Reproduced with permission from Ha et al., Small 19, 2300240 (2023). Copyright 2023 Authors, licensed under a Creative Commons Attribution (CC BY) license.

D. Perovskite solar cells

For perovskite solar cell systems, a number of theoretical studies associated with ML and AIMD simulations have been reported.^189,190 The links connecting the perovskite and layers responsible for carrying electrical charge contain significantly more imperfections compared to those found within the perovskite layer itself—roughly 100 times greater. These imperfections, especially the deep-level ones, noticeably reduce the effectiveness of these devices in converting power.^{72,73,191–193} Recent efforts to reduce these flaws at the connections have primarily focused on enhancing the surface.¹⁹⁴ However, improving the surface of the perovskite that connects with the layer responsible for transporting electrons presents challenges, as the agents used for treating the electron-transporting layer might dissolve during application to the perovskite thin film. Conversely, if a seamless connection between the electron-transporting and perovskite layers could be established, concerns about flaws at the connection might become less significant.

The development of an intermediary layer between a SnO₂ layer that transports electrons and a halide perovskite layer that absorbs light involves combining Cl-bonded SnO₂ with a Cl⁻ containing precursor for the perovskite. This intermediary layer demonstrates coherence at the atomic level, improving the extraction and movement of charges from the perovskite layer while exhibiting fewer imperfections at the connection. This coherent intermediary layer has facilitated the creation of perovskite solar cells achieving an efficiency of 25.8% under standard light exposure. Furthermore, even after continuous exposure to light for 500 h, uncovered devices maintained $\approx 90 %$ of their initial efficiency. These findings provide valuable insights into crafting interfaces that minimize imperfections between metal halide perovskites and layers responsible for transporting electrons.

The experimental analysis suggested the ability to form stable FA(MA)SnCl_x structures between the FAPbI₃ perovskite and SnO₂ without significant distortion in the lattice. However, when TiO₂ was utilized as an electrode, establishing a stable interface proved challenging due to the substantial difference in electronegativity between Ti and Sn [Fig. 26(b)]. As per DFT simulations, the interface between the perovskite and TiO₂ did not demonstrate stability in the TiO₂/FAPbI₃ scenario [Fig. 26(c)]. Conversely, the perovskite and SnO₂ interface exhibited the creation of FASnCl₃ interlayers at the SnO₂/FAPbI₃ interface [Fig. 26(d)], resulting in a favorable alignment of energy bands between SnO₂ and FAPbI₃ [Fig. 26(a)]. The conduction band minimum (CBM) of SnO₂ was slightly lower (by 0.21 eV) than that of FAPbI₃ due to the presence of short strong hydrogen bonding (SSHB) between interfacial oxygen and the FA molecule, stabilizing the dangling bond of the interfacial oxygen.¹⁹⁵

FIG. 26.

View large Download slide

Interlayer formation from Cl-BSO and Cl-CPP. (a) (left) Layered stack of the ETL, perovskite, and HTL with an interlayer between the perovskite absorber and ETL, and (right) the energy diagrams of the ETL, perovskite, and HTL. (1): carrier transport to the contact interfaces via the interlayer. (2): back recombination of the extracted carriers, which remained in the interface region, will be inhibited due to the coherent interlayer. (b) Analysis of the Cl ions contained in the SnO₂ and TiO₂ electrodes. (c) Initial SGPR-based ML DFT and then refined DFT simulation of the formation of the interlayer between the perovskite and SnO₂. (d) and (e) DFT optimized geometries of (d) SnO₂/FAPbI₃ with stable interfacial FASnCl₃ perovskite monolayer and (e) TiO₂/FAPbI₃ without stable interface. [Pb (black), I (purple), Cl (green), C (brown), N (light blue), H (white), Sn (dark blue), Ti (blue), and O (red)]. Figure adapted from Ref. 191. Reproduced with permission from Min et al., Nature 598, 444 (2021). Copyright 2021 Authors, under exclusive license to Springer Nature Limited.

E. Universal first-principles force-fields for hydrocarbon systems

Section V E presents the development of universal ML IAPs for a wide range of hydrocarbons, encompassing alkanes, alkenes, and aromatic compounds. Leveraging the SGPR algorithm and on-the-fly adaptive sampling, these potentials were generated based on the PBE+D3 level of DFT. The ML models were trained on diverse hydrocarbon molecules and their clusters encompassing intra- and inter-molecular interactions. Notably, the combined potentials accurately capture the intricate potential energy surfaces for various hydrocarbon systems, demonstrating excellent transferability from small molecules to long alkanes and polyenes and complex ring-chain molecular systems. Furthermore, the ML-IAPs successfully predicted the structures and energies of β-carotene monomer and dimer. Simulation outcomes for liquid ethylene and toluene crystals showcased the ML IAPs' ability to reproduce molecular volumes and crystalline stability differences between different phases. These findings signify the potential of these ab initio-level force fields to advance toward modeling universal organic, polymeric, and biomolecular systems. Classical force-fields like OPLS,^196–200 AMBER,²⁰⁷ GROMACS,²⁰² and CHARMM²⁰³ have long served as essential tools for simulating organic, biological, polymeric, and macro-molecular systems. Though these force-fields reproduce some experimental data such as conformation reasonably well, they fall short in describing non-covalent inter-molecular interactions^194–199 as well as helicity.^210–214 It is particularly difficult to accurately describe both structures and properties.^215,216 Furthermore, their reliance on descriptors such as bond stretching and angle bending limits their applicability to non-reactive systems, hindering their ability to capture diverse molecular behaviors accurately.

In contrast, ML IAPs offer non-parametric representations, but their efficacy depends on comprehensive descriptors that account for many-body correlations while remaining invariant to system transformations. ML-IAPs rely on “descriptors” derived from atomic positions for their expression. While bonds, bond angles, and torsional angles can serve as useful descriptors, their utility is limited to non-reactive systems. Consequently, there is a necessity for more versatile descriptors that capture many-body correlations, remaining unchanged under the physical system's translations, rotations, and permutations of identical atoms.¹²⁸ Ideally, a seamless correspondence should exist from the coordinate space to the descriptor space. If the transformation of descriptors lacks injectivity, it can impose a lower boundary on the achievable accuracy of ML-IAPs.²¹⁷ Thus, the ongoing challenge involves developing more comprehensive descriptors.

Sampling data for ML-IAPs pose a significant challenge, necessitating training data from various potential energy surface domains. Most of these ML-IAPs are often generated on-the-fly with MD simulations.^{43,44,54,91,101} Active or on-the-fly learning requires assessing the model's uncertainty at given input points. Bayesian inference methods straightforwardly predict uncertainty, while NNs emulate it as variance among predictions from a committee of NNs trained with different architectures or random seeds. However, generating the initial NNs demands substantial effort. While NNs are advantageous when abundant ab initio data are available, the challenge remains in sampling these data across all relevant domains of the potential energy surface. Conversely, kernel-based or Bayesian methods enable on-the-fly model generation. The difficulties in utilizing NNs for on-the-fly training stem from their lack of uncertainty prediction comparable to Bayesian methods and their extensive hyper-parameterization. Re-optimizing all hyper-parameters upon sampling new data becomes cumbersome in high-dimensional spaces and is less conducive to on-the-fly training. Among Bayesian methods, exact GPR offers optimal accuracy but scales computationally as $O (n^{3})$ with the training data's size (⁠ $= n$ ⁠). In contrast, SGPR scales more favorably as $O (n)$ while retaining the favorable characteristics of GPR. Hence, SGPR emerges as an ideal choice for on-the-fly sampling of optimal datasets in generating ML IAPs.

The focus of this study lies in developing ML IAPs using SGPR for a broad spectrum of hydrocarbons, encompassing alkanes and alkenes with different structures and different degrees of unsaturation. Though triple bonds have not been included yet, here the work has considered not only saturated, olefinic, and aromatic hydrocarbons, which include alkane, alkene, alkadiene, and alkatriene molecules (having zero to arbitary number of C=C double bonds) with branches but also cyclic saturated/unsaturated and aromatic hydrocarbons, and further their combined structures.^81,82 This research aims to lay the groundwork for a universal ab initio force-field specific to this subset of hydrocarbons. The goal is to bridge the accuracy gap observed in conventional force-fields and create a more versatile model capable of accurately describing diverse organic, bioorganic systems as well as polymers and macromolecular systems.

1. Machine learning potentials

ML IAPs adopt a holistic approach by encapsulating all atoms within a specified cutoff radius as descriptors. For any configuration of N atoms represented as x, a compilation of descriptors

x = {ρ_{i}}_{i = 1}^{N}

is created. Each ρ_i relies on the relative coordinates of all atoms in the LCE of atom i within the cutoff radius. The potential energy, E(x), is then computed:

E (x) = \sum_{i = 1}^{N} E (ρ_{i}),

(147)

where

E

is a fictional energy representing the interaction of atom i with its neighboring atoms. Kernel-based regression methods define

E

as

E (ρ) = \sum_{j = 1}^{m} K (ρ, χ_{j}) w_{j} = k_{ρ m} w,

(148)

where the matrices are expressed with bold fonts,

z = {χ_{j}}_{j = 1}^{m}

is the set of reference/inducing descriptors, w is the vector of weights for the inducing descriptors, and

K

is a similarity/covariance kernel. Namely, the potential energy E(x) is a sum over atoms i, involving a similarity/covariance kernel

K

operating on the descriptors ρ_i and inducing descriptors χ_j, weighted by w_j. These weights are determined to accurately replicate the potential energy and forces of an ab initio dataset

X = {x_{k}}_{k = 1}^{n}

⁠. In SGPR,⁵⁴^,

w = {(σ^{2} k_{m m} + k_{n m}^{T} k_{n m})}^{- 1} k_{n m}^{T} Y,

(149)

where

k_{m m}

and

k_{n m}

are the inter-inducing and data-inducing covariance matrices, respectively, Y are the data potential energies (and forces), and σ is the noise hyper-parameter. The noise scale σ and other possible hyper-parameters in the kernel

K

are optimized to maximize likelihood of the data energies. For kernel similarity, a variation of SOAP¹²⁸ defined in Ref. 54 is employed. The inducing descriptors are selected from existing data or sampled during on-the-fly learning in MD simulations, aiding in uncertainty prediction and facilitating the learning process. With Bayesian methods, the predictive variance serves as a threshold for sampling new data or inducing descriptors. In SGPR, the predictive variance for

E (ρ)

is given by

ν (ρ) = k_{ρ ρ} - k_{ρ m} k_{m m}^{- 1} k_{ρ m}^{T} + σ^{2} k_{ρ m} {(σ^{2} k_{m m} + k_{n m}^{T} k_{n m})}^{- 1} k_{ρ m}^{T} .

(150)

For calculation of w from Eq. , the matrix inversion

{(σ^{2} k_{m m} + k_{n m}^{T} k_{n m})}^{- 1}

can be bypassed by converting this equation into a linear system.^54,217 However, in Eq. , this matrix inversion can be computationally burdensome²¹⁸ for on-the-fly learning, prompting the use of approximations for computational effectiveness in spite of occasionally slightly sub-optimal sampling outcomes.

In this work, for all first-principles calculations, the Vienna ab initio simulation package (VASP)¹⁰⁷ using the Perdew–Burke–Ernzerhof (PBE) generalized-gradient approximation (GGA) functionals¹⁰⁹ with correction method of Grimme et al.²¹⁹ for van der Waals interactions (D3) was used. The python package $AUTOFORCE$ ¹⁰⁵ for on-the-fly generation of SGPR models.⁵⁴ was used. For the accuracy of ML IAPs, the coefficient of determination denoted by $R^{2} = 1 - \sum_{i} {(f_{i} - {\tilde{f}}_{i})}^{2} / \sum_{i} {(f_{i} - \bar{f})}^{2}$ is reported, where ${f_{i}}$ are the first-principles (PBE+D3) forces, $\bar{f}$ is their average, and ${{\tilde{f}}_{i}}$ are the ML forces.

Here, the training with a particular molecule means MD simulations with an adaptive sampling algorithm similar to Ref. 54 for on-the-fly ML of the potential energy surface. Accurate variance prediction is central for efficient on-the-fly training. All available conformers of a given molecule are used for training. Each conformer is simulated for 3–6 ps of MD at 300 K. Molecules are queued for training based on their simplicity. Owing to the abundance of relevant molecules, trained molecules are split into several groups and independent SGPR models are trained in parallel using the AutoForce package. In this study, we considered expert models for linear or branched systems (C1–C8 alkane, C3–C10 isoalkane, C2–C6 alkene, C3–C6 branched alkene, C3–C8 alkadiene, C7–C9 alkatriene), cyclic systems (C3–C10 cycloalkane, C3–C8 cycloalkene, C4–C10 methyl-/dimethyl-cycloalkane, bicycloalkane), and aromatic systems (benzene, toluene, xylene, trimethylbenzene, mesitylene) to build the universal model. To consider the intermolecular interactions, expert models for liquid phases of alkane (methane, ethane, propane, butane), alkene (ethylene), and aromatic (benzene) are trained also and denoted as alkane-inter, alkene-inter, and aromatic-inter, respectively. The model for alkane is extended by simulations for small clusters and liquid phases of methane, ethane, propane, and butane for learning inter-molecule interactions. This extended alkane model will be distinguished as “alkane^†.” Later, independent MD simulations are repeated for sampling independent testing sets for all conformers of the above molecules.

2. Macromolecular systems

a. Expert models

The expert SGPR model for each category of molecules is generated on-the-fly by short MD simulations (3–6 ps) of the associated molecules at 300 K. As an example, training for the elementary alkane molecules is shown in Fig. 27. The model trained after each MD simulation is used as the initial model for the next. In this way, the forces for larger molecules are inferred from the data for smaller ones more frequently and less ab initio calculations are performed for on-the-fly training. Expert models are also trained for cyclic hydrocarbons (cycloalkane, cycloalkene, and bicycloalkane) and aromatic systems, and then are combined with the previously generated model for alkane, isoalkane, alkene, and alkadiene. Each generated expert model is validated by comparing it with the DFT calculation results. Similar MD simulations are repeated without on-the-fly training for generating independent testing sets. The performance of the expert models on their testing sets are shown in Table III by means of mean absolute error (MAE) and R². These SGPR models closely approximate the ab initio PBE+D3 potential energy surface in their respective domains.

FIG. 27.

View large Download slide

Initial steps for on-the-fly ML MD for elementary alkane molecules. Red bullets indicate the ab initio potential energy of the configurations sampled by the SGPR model. Green bullets are for single-point ab initio tests.⁸¹ Figure taken from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

TABLE III.

Testing the expert models in their domains. N_test is the number of test samples. Energy^* is the potential energy per-atom.⁸¹ Adapted from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

Group	N_test	Energy^* MAE	Force MAE	Force R²
		(meV)	(eV/Å)	Testing/training
Alkane	448	4.9	0.076	0.978/0.993
Isoalkane	308	3.0	0.080	0.977/0.989
Alkene-1	142	8.1	0.110	0.966/0.985
Alkene-2	88	5.3	0.110	0.955/0.98
Alkene-3	238	4.0	0.087	0.975/0.991
Alkadiene	374	4.5	0.110	0.959/0.978
Alkatriene	238	3.2	0.108	0.992/0.966
Cycloalkane	253	2.3	0.074	0.983/0.990
Cycloalkene	77	3.4	0.110	0.963/0.984
Methyl/dimethyl-cycloalkane	946	5.2	0.092	0.970/0.989
Bicycloalkane	594	2.1	0.075	0.980/0.987
Aromatic	83	2.1	0.074	0.997/0.983

Group	N_test	Energy^* MAE	Force MAE	Force R²
		(meV)	(eV/Å)	Testing/training
Alkane	448	4.9	0.076	0.978/0.993
Isoalkane	308	3.0	0.080	0.977/0.989
Alkene-1	142	8.1	0.110	0.966/0.985
Alkene-2	88	5.3	0.110	0.955/0.98
Alkene-3	238	4.0	0.087	0.975/0.991
Alkadiene	374	4.5	0.110	0.959/0.978
Alkatriene	238	3.2	0.108	0.992/0.966
Cycloalkane	253	2.3	0.074	0.983/0.990
Cycloalkene	77	3.4	0.110	0.963/0.984
Methyl/dimethyl-cycloalkane	946	5.2	0.092	0.970/0.989
Bicycloalkane	594	2.1	0.075	0.980/0.987
Aromatic	83	2.1	0.074	0.997/0.983

b. Long alkane molecules

Using the alkane SGPR model, which is trained using only short molecules, fails to reproduce the potential energy surface for longer molecules when they fold and form a coil. When a long molecule folds, different parts of the same molecule interact with each other, which is not reproducible with short molecules. Instead, alkane^† model, which is extended for inter-molecule interactions using small alkane molecules, is applicable. Here we used this model for MD simulations (without on-the-fly training) of long n-alkane molecules for 100 ps at 300 K. Snapshots of C₆₄H₁₃₀ are shown in Fig. 28. A testing set is generated by collecting snapshots every 1 ps. Performance of the alkane^† model is summarized in Table IV, which shows a high degree of transferability. This table also reconfirms that the computational cost is $\sim N^{3}$ for DFT but $\sim N$ for ML.

FIG. 28.

Folding of C64H130 at 300 K. Snapshots are shown every 2 ps.81 Figure taken from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

View large Download slide

Folding of C₆₄H₁₃₀ at 300 K. Snapshots are shown every 2 ps.⁸¹ Figure taken from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

TABLE IV.

Testing alkane^† SGPR model on long alkane molecules. $T$ indicates the elapsed time (seconds) for a single-point energy and force calculation on a node with 20 CPUs. Energy^* is the potential energy per-atom.⁸¹ Adapted from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

Molecule	$T_{M L}$	$T_{DFT}$	Energy^* MAE	Force MAE	R²
	(s)	(s)	(meV)	(eV/Å)	(testing)
C₈H₁₈	0.44	82.07	3.3	0.070	0.985
C₁₆H₃₄	0.65	282.36	4.7	0.085	0.977
C₃₂H₆₆	1.08	1836.88	4.2	0.091	0.973
C₆₄H₁₃₀	1.92	12310.9	4.8	0.095	0.970

Molecule	$T_{M L}$	$T_{DFT}$	Energy^* MAE	Force MAE	R²
	(s)	(s)	(meV)	(eV/Å)	(testing)
C₈H₁₈	0.44	82.07	3.3	0.070	0.985
C₁₆H₃₄	0.65	282.36	4.7	0.085	0.977
C₃₂H₆₆	1.08	1836.88	4.2	0.091	0.973
C₆₄H₁₃₀	1.92	12310.9	4.8	0.095	0.970

c. Universal model

The SGPR models for alkane, isoalkane, alkene-(1, 2, 3), and alkadiene are combined into a single universal model for linear/branched hydrocarbons with arbitrary number of double bonds. The model was rebuilt using all of the ab initio data and inducing descriptors, which were sampled during training of the expert models.⁶³ As shown in Table V, the universal model is more accurate than the expert models even in their own domains. The correlation of the predicted energies with the exact ab initio energies for all testing sets is shown in Fig. 29.

TABLE V.

Testing the universal model in all domains. The testing sets are identical to those in Table III.⁸¹ Adapted from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

Group	Energy^* MAE	Force MAE	Force R²
	(meV)	(eV/Å)
Alkane	3.4	0.055	0.989
Isoalkane	3.3	0.060	0.986
Alkene-1	5.6	0.075	0.98
Alkene-2	5.1	0.079	0.973
Alkene-3	4.1	0.071	0.983
Alkadiene	4.4	0.081	0.977

Group	Energy^* MAE	Force MAE	Force R²
	(meV)	(eV/Å)
Alkane	3.4	0.055	0.989
Isoalkane	3.3	0.060	0.986
Alkene-1	5.6	0.075	0.98
Alkene-2	5.1	0.079	0.973
Alkene-3	4.1	0.071	0.983
Alkadiene	4.4	0.081	0.977

FIG. 29.

Predicted energies of the universal model for all testing sets.81 Figure taken from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

View large Download slide

Predicted energies of the universal model for all testing sets.⁸¹ Figure taken from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

For a case study, several linear and branched isomers of C₁₁H₁₈ were considered. There are 3C=C double bonds in these isomers (see Fig. 30). It is worth noting that the training molecules contained at most 2C=C double bonds (alkadiene). These structures are relaxed using the universal ML-IAP model, PBE+D3, B3LYP+D3, and OPLS-AA force-field by minimizing the forces below 0.01 eV/Å. B3LYP+D3 is included since it is a popular functional for molecular systems. Additionally, we are intrigued to check how the error of ML-IAP compares to the disagreement of two widely used density functionals. The energies for the optimized structures, depending on the method used for relaxation, are shown in Fig. 31. For structures optimized using PBE+D3, single-point ML-IAP and OPLS-AA energies are compared. The energies for optimized structures using ML-IAP, PBE+D3, and B3LYP are in reasonable agreement while OPLS-AA fails to reproduce the trend. It should be noted that OPLS-AA is an empirical potential designed to statistically reproduce some experimental observations rather than single-point relaxed structures.

FIG. 30.

Isomer structures of C11H18 for testing the universal model. 15 different isomers are prepared for C11H18.81 Figure taken from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

View large Download slide

Isomer structures of C₁₁H₁₈ for testing the universal model. 15 different isomers are prepared for C₁₁H₁₈.⁸¹ Figure taken from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

FIG. 31.

Relative energies ΔE=Ei−E1 (i = 1–15, see Fig. 30) of C11H18 isomers relaxed with PBE+D3, ML-IAP, B3LYP+D3, and OPLS-AA. (sp) are single-point energies for configurations relaxed by PBE+D3.81 Figure taken from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

View large Download slide

Relative energies $Δ E = E_{i} - E_{1}$ (i = 1–15, see Fig. 30) of C₁₁H₁₈ isomers relaxed with PBE+D3, ML-IAP, B3LYP+D3, and OPLS-AA. (sp) are single-point energies for configurations relaxed by PBE+D3.⁸¹ Figure taken from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

d. Rotation of CH₃ group in ethane

The accuracy of universal ML-IAPs generated by combining expert models is tested first. As a simplest case, nudged elastic band (NEB) calculation is conducted to predict the rotation barrier of CH₃ group in ethane. The barrier to rotate CH₃ group for 120° is calculated using universal ML-IAPs and compared with PBE+D3 results. As a result, the barrier of the reaction was predicted to be 0.12 eV, which is almost the same with the PBE+D3 results showing 0.11 eV.

e. Dimerization of β-carotene

The universal ML-IAPs are used to predict the structures of molecular systems containing various hydrocarbon groups. Here, the structures of β-carotene monomer and dimer are investigated using the universal potential. β-carotene is an organic molecule having eight isoprenes and two trimethyl-cyclohexenes in its structure. The expert models including C–C single bonds, C=C double bonds, and cyclic hydrocarbons are used to construct the universal model. The relative energies of 20 different conformers of β-carotene monomer are predicted to be close to the PBE+D3 results within a low RMSE of 0.12 eV.

f. Liquid ethane and ethylene

The alkane^† model is used for simulations of liquid ethane at the experimental boiling temperature. In all, 216 molecules are used for MD simulations in NPT ensemble with constant pressure of 1 bar. The duration of simulation was 100 ps for equilibration and 200 ps for averaging. Heat of vaporization is calculated as

Δ H = H_{g} - H_{l} \approx (u_{g} + R T) - (u_{l} + p v_{l}),

(151)

where u_g/v_g and u_l/v_l are the internal energy/volume per molecule in the gas and liquid phases, respectively, p is the pressure, and H_g is calculated with ideal gas approximation (pv_g = RT). The results are compared with the experimental data in Table VI. In this case, although ML-IAP yields v_l in good agreement with experiment, it slightly overestimates

Δ H

⁠. However, unlike empirical force-fields, this ML-IAP is truly a first-principles potential whose accuracy is limited by PBE+D3 level of DFT. Note that in this case a direct testing of DFT, if not impossible, would be extremely difficult. Thus, surrogate ML-IAP models can be used for testing expensive, lower level, first-principles potentials against experimental observations.

TABLE VI.

Simulation of liquid ethane at 184.52 K and 1 bar. v_l is the volume per molecule in the liquid phase. $Δ H$ is the heat of vaporization. Adapted from Ref. 81. Reproduced with permission from Hajibabaei et al., J. Phys. Chem. A 125, 9414 (2021). Copyright 2021 Authors, licensed under a Creative Commons Attribution (CC BY) license.

	Units	ML (alkane^†)	Experimental
v_l	(Å³)	90.3(4)	91.5
$Δ H$	(kcal/mol)	4.01(2)	3.52

Also, the molecular volume of liquid ethylene is predicted using ML-IAPs. The NPT MD simulations of 34 ethylene molecules are conducted for 30 ps with a constant pressure 1 bar and temperature 169.5 K. The volume per molecule v_l is predicted to be 76 ± 6 Å³ in good agreement with the experimental value²²⁰ of 81.5 Å³. Some discrepancy from the experimental value seems to arise from the inherent discrepancy in DFT-D3 calculation.

g. Toluene crystal

The prediction of crystal structures is accelerated by using the ML potential trained with DFT potential energies/forces. Here, the universal model predicts the stable toluene crystal structure (Fig. 32). The α-phase is predicted to be the most stable state in crystal, while the β-phase is metastable state.²²¹

FIG. 32.

Crystal structures of (a) α-toluene and (b) β-toluene.82 Figure taken from Ref. 82. Reproduced with permission from Hajibabaei et al., ACS Phys. Chem. Au 2, 260–264 (2022). Copyright 2022 Authors, licensed under a Creative Commons Attribution (CC BY) license.

View large Download slide

Crystal structures of (a) α-toluene and (b) β-toluene.⁸² Figure taken from Ref. 82. Reproduced with permission from Hajibabaei et al., ACS Phys. Chem. Au 2, 260–264 (2022). Copyright 2022 Authors, licensed under a Creative Commons Attribution (CC BY) license.

3. Discussion

This study outlines the comprehensive generation of a universal ML-IAP designed to predict molecular and crystal structures of hydrocarbons. This model encapsulates various hydrocarbon types, including alkanes, unsaturated alkenes with diverse C=C double bond numbers, aromatic compounds, open-chain, and cyclic compounds. Using the SGPR algorithm and on-the-fly learning, the model emerges from less than 500 atomic configurations of elementary molecules, along with clusters and liquid phases of small molecules to incorporate inter-molecule interactions.

Expert models tailored to distinct hydrocarbon groups are initially developed through on-the-fly ML and then consolidated into the universal ML-IAPs. These universal models, originating from the expert models trained for 90 molecules (350 conformers) and five bulk phases, exhibit exceptional transferability across varied hydrocarbon systems. They effectively handle complex scenarios like the folding and coiling of long alkane molecules, showcasing reliability even under such structural changes. Simulations of liquid ethane at its boiling temperature yield results reasonably aligned with experimental observations.

The model's prowess extends to predicting the minimum energy structure of C₁₁H₁₈, displaying impressive agreement with direct ab initio simulations while outperforming traditional force-fields like OPLS-AA. Although OPLS-AA is empirically successful in replicating certain experimental observations, ML-IAPs, being surrogates for expensive first-principles potentials, provide more adaptable avenues for studying bond breaking/forming/rearrangements reaction-paths involving atom transfer, transition-states, and phase transitions. The universal ML-IAPs accurately predict the structure and interaction energies of complex compounds like the $β$ -carotene dimer and effectively describe molecular volume and phase transition, as exemplified in the investigation of the toluene crystal's phases. The significance of this work lies in efficiently generating a universal model from separately trained atomic potentials for various subgroups, enabling its broad applicability across multi-component systems—a feat unattainable for other ML methods. It showcases one of the most extensive explorations of hydrocarbon phase spaces, employing fewer than 1000 configurations as the basis for SGPR potential. Employing a Bayesian criterion during on-the-fly MD simulations allows automatic discarding of configurations with redundant information. It emphasizes that data quality, in terms of diversity and completeness, significantly impacts the accuracy and transferability of ML potentials, surpassing the importance of data quantity. Moreover, it explores the potential of using more comprehensive descriptor schemes for improved accuracy and introduces intriguing concepts like “recursively embedded atom NNs” for incorporating higher-order correlations. It also highlights the importance of balancing global exploration and local exploitation in data sampling for developing accurate ML-IAPs.

The SGPR algorithm emerges as a robust and rapid tool for on-the-fly ML of potential energy surfaces, offering excellent transferability despite being based on a relatively small number of configurations. It introduces a modular approach through expert models, allowing the combination of diverse models into a singular global model. The study underscores the significance of building efficient yet globally diverse datasets, a crucial step in the automation and optimization of ML-IAPs.

VI. PERSPECTIVES AND CONCLUDING REMARKS

ML has gained recognition for its precision in modeling quantum potential surfaces and forces. While NNs are frequently employed for this purpose, they require frequent updates from numerous trajectories. GPR provides an alternative. However, traditional GPR models still necessitate large datasets for training due to extensive covariance matrix elements. This complicates continuous learning, particularly during phase transitions and chemical reactions. SGPR addresses this issue by significantly reducing the number of parameters, which allows for smaller training datasets. SGPR's ability to perform rank reduction improves the computational performance by several orders of magnitude (see Table VII). Moreover, by utilizing path-independent feature parameters through active learning instead of traditional atomic coordinates, the SGPR-ML based algorithm facilitates the creation of universal quantum mechanical potentials for complex systems, establishing a foundational platform for virtual material modeling and broadening applicability across various quantum mechanical systems. Path-independent parameters can be compiled into a built-in library, enabling their reuse in related systems without the need for re-training. The transferability of SGPR's covariance matrix means it remains effective without updates. This facilitates the creation of a comprehensive covariance matrix library through cumulative learning, enhancing the efficiency and speed of simulations for complex material systems. This makes SGPR versatile for the simulations of a range of small to complex quantum mechanical systems at the levels of advanced theories like coupled cluster theory with perturbative triple excitations [CCSD(T)],²²² matrix product states (MPS),^223,224 and quantum Monte Carlo (QMC) simulations^225,226 alongside DFT.

TABLE VII.

Computational cost of normal GPR vs SGPR. n is the number of snapshots in the data, N is the number of atoms in each snapshot, n_c is the average number of atoms inside a sphere with cutoff radius r_c, and m is the number of inducing LCEs sampled by SGPR. In the last column, an exemplary ratio in computational cost for GPR over SGPR is estimated for the case where n = N = 100, n_c = 33, and m = 1000, similar to the cases that we have often encountered.

Method	DFT-data	Accuracy	Paper
Cost	GPR	SGPR	GPR/SGPR
Matrix inversion	(1 + 3N) $^{3} n^{3}$	(1 + 3N)nm²	∼10³
Kernel construction	(1 + 3n_c) $^{2} N^{2} n^{2}$	(1 + 3n_c)Nnm + m²	∼10³

Method	DFT-data	Accuracy	Paper
Cost	GPR	SGPR	GPR/SGPR
Matrix inversion	(1 + 3N) $^{3} n^{3}$	(1 + 3N)nm²	∼10³
Kernel construction	(1 + 3n_c) $^{2} N^{2} n^{2}$	(1 + 3n_c)Nnm + m²	∼10³

SGPR's rank reduction, empowered by selecting LCE features via similarity measures in active learning, replaces atomic coordinates with characteristic geometric features. This method helps achieve chemical accuracy across various spatial configurations of quantum mechanics potentials. By aggregating covariance elements from subsystems and analyzing their interrelationships, SGPR is capable of generating universal quantum mechanical potentials, enabling integrated first-principles calculations and ML simulations on a comprehensive platform for virtual materials reality. This advanced capability is generally unachievable with most traditional ML methods. Meanwhile, the covariance matrices created through SGPR's rank reduction are transferable, meaning that their updates are not required even during events such as phase transitions, molecular transformations, or chemical bond formation/breaking. This feature allows the combination of systems, where inter-system covariance matrices are integrated gradually, building a comprehensive library. This library then facilitates rapid and precise simulations for complex systems by simply retrieving relevant datasets, effectively realizing the structural and phenomenological characteristics of these systems. SGPR thus distinguishes itself with several key features: rapid training due to minimal LCE-based rank reduction, high-speed computations suitable for large or complex systems, and the generation of a universal, transferable potential that supports both expansive learning across complex systems and the development toward self-learning super-AI or digital twin. By merging SGPR with active learning, meta-dynamics, and quantum mechanical calculations, it would be possible to develop a self-improving AI and digital twin capable of innovating new materials and devices. SGPR offers a reduction in computational complexity and biases while ensuring quicker, more accurate simulations for sizable and intricate systems. Its unique attributes, such as linear scalability and adaptability to varied systems, represent a substantial leap over conventional ML methodologies, encouraging creative and effective approaches in materials science simulations and designs. The strengths of SGPR include (i) fast training due to reduced parameter needs; (ii) high-speed computation suitable for large and complex systems (linear scaling); (iii) chemical accuracy and the ability to produce universal potentials; and (iv) scalability through cumulative learning, forming a base for a self-learning AI or Digital Twin.

This approach has been demonstrated in the simulation of superionic materials for fast-charging batteries and electro-catalytic reactions, showing significant improvements over traditional methods. SGPR's low-rank approximation of the covariance matrix boosts computational performance dramatically, which we have quantified in comparative studies. Additionally, SGPR integrates seamlessly with on-the-fly ML methods, eliminating the need for constant retraining of ML models during simulations. This makes it particularly suited for extensive materials studies involving solid electrolytes, catalysts, solar cells, and even complex biological molecules, where it enables rapid simulation and exploration of potential materials without the computational overhead typically associated with such tasks. Thus, SGPR represents a transformative step forward in the use of ML for quantum mechanical simulations, providing a robust platform for the development of new materials and the detailed study of existing ones under various conditions, all while reducing the computational expense and increasing the scope of potential applications.

In this work, we reviewed the applications of SGPR for extensive materials involving bulk materials, solid electrolytes, batteries, catalysts, solar cells, and macromolecular systems. Leveraging this tool allows to bypass traditional computational constraints, enabling extensive exploration of material compositions. This approach is expected to significantly contribute to the development of safer and more efficient solid electrolytes for future energy storage applications. Here, as an exemplary application, we outline a technical strategy and the expected outcome of SGPR as a comprehensive virtual lab capable of simulating and predicting properties of electrolytes across various systems, significantly advancing our capability to design efficient, stable, and safe solid electrolytes for future energy storage solutions.

The Innovative Approach will push the boundaries of materials science by integrating advanced ML techniques with traditional chemical physics approaches to develop more advanced materials, which is crucial for designing next-generation novel materials.

SGPR-based ab initio molecular dynamics (AIMD): AIMD excels in simulating complicated systems. However, AIMD's computational demands scale cubically with the number of atoms, making it feasible only for smaller systems and shorter durations. SGPR-enhanced AIMD can handle larger atom counts over extended periods, demonstrating its potential for larger, more complex simulations.
Integration of quantum mechanics and ML: This SGPR-based approach ushers in a new era for modeling with significant applications such as:
- Learning from existing data to develop strategies, similar to how AlphaGo learned to master Go toward AlphaGo-0.
- Generating substantial data on unknown material properties using quantum mechanical computations.
- Creating a universal quantum mechanics potential that combines quantum mechanics and ML, enabling rapid computations.
- Exploring structural space and phase transitions through meta-dynamics.
- Developing self-learning systems that evolve into Super AI and Digital Twins, capable of innovating new materials and devices.
Challenges with conventional ML methods: Traditional GPR methods require inverting large covariance matrices, which is impractical for systems with thousands of trajectories due to computational constraints. SGPR addresses this by reducing feature dimensions, which simplifies the learning process and enhances the speed and scalability of simulations.
Key features of SGPR:
- Fast training and superfast computing capabilities suitable for extensive systems due to linear scaling.
- Ability to generate transferable universal potentials that facilitate rapid, large-scale quantum calculations.
- Minimal data requirements for accurate potential generation, leveraging the cumulative learning of key features.

To conclude, by merging computational power with innovative ML strategies, this research stands to significantly impact the field of energy storage, contributing to safer and more reliable battery technologies and other novel functional materials including catalysts and solar cells that are essential for a range of modern applications.

ACKNOWLEDGMENTS

This research was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT: Ministry of Science and ICT) (No. RS-2024-00346451) and the Program for Cultivating the Next Generation of Outstanding Scientists (No. 1.240033.01) by UNIST Basic Science Institute.

AUTHOR DECLARATIONS

Conflict of Interest

The authors have no conflicts to declare.

Author Contributions

Soohaeng Yoo Willow, Amir Hajibabaei, Miran Ha, David ChangMo Yang, Chang Woo Myung, Seung Kyu Min, and Geunsik Lee contributed equally to this work.

Soohaeng Yoo Willow: Writing – original draft. Amir Hajibabaei: Writing – original draft. Miran Ha: Writing – original draft. David ChangMo Yang: Writing – review & editing. Chang Woo Myung: Writing – review & editing. Seung Kyu Min: Writing – review & editing. Geunsik Lee: Writing – review & editing. Kwang S. Kim: Conceptualization; Funding acquisition; Project administration; Writing – original draft; Writing – review & editing.

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

APPENDIX: MATHEMATICAL BACKGROUND

1. Mathematical expressions

a. Matrix algebra

Matrix inversion lemma:

The matrix inversion lemma shows

\begin{matrix} {(A + x)}^{- 1} & = {[A (1 + A^{- 1} x)]}^{- 1} = {(1 + A^{- 1} x)}^{- 1} {(1 + A^{- 1} x) - A^{- 1} x} A^{- 1} \\ = A^{- 1} - {(1 + A^{- 1} x)}^{- 1} A^{- 1} x A^{- 1} = A^{- 1} - {[x^{- 1} A (1 + A^{- 1} x)]}^{- 1} A^{- 1} \\ = A^{- 1} - {[(x^{- 1} + A^{- 1}) A]}^{- 1} A^{- 1} = A^{- 1} - A^{- 1} {(A^{- 1} + x^{- 1})}^{- 1} A^{- 1} . \end{matrix}

(A1)

Therefore,

{(A + B)}^{- 1} = A^{- 1} - A^{- 1} {(A^{- 1} + B^{- 1})}^{- 1} A^{- 1}

or

{(A^{- 1} + B^{- 1})}^{- 1} = A - A {(A + B)}^{- 1} A = B - B {(A + B)}^{- 1} B .

Consequently, the Woodbury, Sherman, and Morrison formula shows

{(A + U C V^{T})}^{- 1} = A^{- 1} - A^{- 1} U {(C^{- 1} + V^{T} A^{- 1} U)}^{- 1} V^{T} A^{- 1} .

(A2)

A useful expression is

{(σ^{2} I + Q Q^{T})}^{- 1} = σ^{- 2} I - σ^{- 2} Q {(σ^{2} I + Q^{T} Q)}^{- 1} Q^{T} .

(A3)

Matrix determinant lemma:

For the dyadic product

a b^{T}

of a column vector

u

and a row vector

v^{T}

⁠, the matrix determinant lemma satisfies

[\begin{matrix} I & 0 \\ v^{T} & 1 \end{matrix}] [\begin{matrix} I + u v^{T} & u \\ 0 & 1 \end{matrix}] [\begin{matrix} I & 0 \\ - v^{T} & 1 \end{matrix}] = [\begin{matrix} I & u \\ 0 & 1 + v^{T} u \end{matrix}] .

(A4)

Therefore,

Sparse Gaussian process based machine learning first principles potentials for materials simulations: Application to batteries, solar cells, catalysts, and macromolecular systems

I. INTRODUCTION

II. THEORETICAL BACKGROUND

A. Introduction to machine learning (ML) methods

B. Notations

C. Representation of the potential energy

D. Complexity of machine learning models

E. Neural networks

F. Gaussian processes

1. Multivariate normal distributions

2. Gaussian process regression

3. Hyper-parameters

4. Covariance matrix for energies and forces

G. Potential energy inference

H. Bayesian linear regression

I. Compressed sensing

J. Sparse Gaussian process regression (SGPR)

K. Variational lower bound for log-likelihood

L. Adaptive sampling algorithm

M. AutoForce: A SGPR-based ML package for atomistic simulations

N. Benchmark test of SGPR-based ML: Silicon

O. Universal potentials

P. Sparse Bayesian committee machine potential

Q. Comparison of GNNs and SGPR

R. Summary of Sec. II

III. BASIC FUNCTIONS

A. General aspects

B. Combining descriptors

C. Rotational symmetry

D. Two-body and many-body

E. Coulomb matrix

F. Behler–Parrinello symmetry functions

G. Smooth overlap of atomic positions

H. Summary of Sec. III

IV. DYNAMIC PROPERTIES

A. Conductivity

B. Ionic conductivity in crystals

C. Diffusion in solids

D. Linear response theory

E. Diffusion constant from molecular dynamics

F. Grain boundaries and long-range diffusivity

G. Diffusivity and activation energy

H. Summary of Sec. IV

V. APPLICATIONS: SGPR-ASSISTED MATERIALS SIMULATIONS

A. Sulfide solid electrolytes

1. Li2S-P2S5 system

2. Training

3. Testing

4. Computational performance

5. α and β phases

6. Glass phase

7. Li–X–Y ternary systems

a. Searching algorithm

b. Universal ML interatomic potentials for surveying ternary solid electrolytes

8. Universal potentials: Application to quaternary systems

9. Summary for SGPR simulations of solid electrolytes

B. Lithium ion batteries

1. Al-doping in Li-ion-battery cathode materials

2. DFT calculations

3. Structure prediction

4. Structures of LRNO and Al-LRNO

5. Stability fading of Al-LRNO due to delithiation at high-voltages

6. Enhanced diffusion by Al-doping

7. Electronic origin of stabilization and Discussion

8. Sodium-substituted disordered rock salt cathode: Reduced energy barrier for oxygen evolution

C. Electrocatalysts

1. High-performing atomic electrocatalysts

2. HER/OER/ORR

3. CO2RR

4. NRR

5. ClER

D. Perovskite solar cells

E. Universal first-principles force-fields for hydrocarbon systems

1. Machine learning potentials

2. Macromolecular systems

a. Expert models

b. Long alkane molecules

c. Universal model

d. Rotation of CH3 group in ethane

e. Dimerization of β-carotene

1. Li₂S-P₂S₅ system

3. CO₂RR

d. Rotation of CH₃ group in ethane