Machine-learning of atomic-scale properties based on physical principles

TL;DR

This paper reviews kernel regression methods for modeling atomic-scale properties, introduces the SOAP representation for atomic environments, and discusses advanced generalizations for improved material property predictions.

Contribution

It provides a detailed account of the SOAP representation and kernel, connecting it to other atomic structure descriptors, and discusses recent generalizations for complex property prediction.

Findings

SOAP arises from smooth atomic densities

Kernel methods can predict energies and forces

Generalizations enable tensorial property prediction

Abstract

We briefly summarize the kernel regression approach, as used recently in materials modelling, to fitting functions, particularly potential energy surfaces, and highlight how the linear algebra framework can be used to both predict and train from linear functionals of the potential energy, such as the total energy and atomic forces. We then give a detailed account of the Smooth Overlap of Atomic Positions (SOAP) representation and kernel, showing how it arises from an abstract representation of smooth atomic densities, and how it is related to several popular density-based representations of atomic structure. We also discuss recent generalisations that allow fine control of correlations between different atomic species, prediction and fitting of tensorial properties, and also how to construct structural kernels---applicable to comparing entire molecules or periodic systems---that go beyond an additive combination of local environments.

Equations

$$f(\mathcal{A})=\sum _{\mathcal{B}\in M}{x}_{\mathcal{B}}K(\mathcal{A},\mathcal{B}),$$

$$\ell (\mathbf{x})=\Vert \mathbf{Kx}-\mathbf{y}{\Vert }^2.$$

$$\mathbf{x}={\mathbf{K}}^{-1}\mathbf{y}.$$

$${\mathbf{x}}_M=({\mathbf{K}}_{MN}{\mathbf{K}}_{NM}{)}^{-1}{\mathbf{K}}_{MN}{\mathbf{y}}_N,$$

$$\mathbf{x}=(\mathbf{K}+\mathbf{\Lambda }{)}^{-1}\mathbf{y}.$$

$${\mathbf{x}}_M=({\mathbf{K}}_{MM}+{\mathbf{K}}_{MN}{\mathbf{\Lambda }}^{-1}{\mathbf{K}}_{NM}{)}^{-1}{\mathbf{K}}_{MN}{\mathbf{\Lambda }}^{-1}{\mathbf{y}}_N,$$

$$\Vert \mathbf{Kx}-\mathbf{y}{\Vert }_{{\mathbf{\Lambda }}^{-1}}^2+\Vert \mathbf{x}{\Vert }_{\mathbf{K}}^2,$$

$$f(\mathcal{C})={\mathbf{K}}_{\mathcal{C}M}{\mathbf{x}}_M.$$

$$d^2(\mathcal{A},\mathcal{B})=K(\mathcal{A},\mathcal{A})+K(\mathcal{B},\mathcal{B})-2K(\mathcal{A},\mathcal{B}).$$

$$\mathbf{D}=US{V}^{\top }.$$

$${\pi }_{\mathcal{A}}=\frac{1}{k}\sum _{\xi =1}^k(u_{\mathcal{A}}^{\xi }{)}^2$$

$$E_{\mathrm{tot}}=E(\mathcal{A})+E(\mathcal{B})={\mathbf{K}}_{\mathcal{A}M}{\mathbf{x}}_M+{\mathbf{K}}_{\mathcal{B}M}{\mathbf{x}}_M.$$

$${}^{\Sigma }{\mathbf{K}}_M=\sum _{\mathcal{A}\in I}{\mathbf{K}}_{\mathcal{A}M},$$

$$E_{\mathrm{tot}}={}^{\Sigma }{\mathbf{K}}_M{\mathbf{x}}_M.$$

$$\Vert {}^{\Sigma }\mathbf{Kx}-{\mathbf{E}}_{\mathrm{tot}}{\Vert }^2,$$

$${\nabla }_{\mathcal{A}}f(\mathcal{A})=\sum _{\mathcal{B}\in M}{x}_{\mathcal{B}}{\overleftarrow{\nabla }}_{\mathcal{A}}K(\mathcal{A},\mathcal{B})={\mathbf{K}}_{\nabla \mathcal{A}M}{\mathbf{x}}_M,$$

$$\Vert {\mathbf{K}}_{\nabla NM}{\mathbf{x}}_M-{\mathbf{y}}_{\nabla N}{\Vert }^2,$$

$$f(\mathcal{A})=\sum _{\mathcal{B}\in M}{\mathbf{x}}_{\nabla \mathcal{B}}\cdot {\overrightarrow{\nabla }}_{\mathcal{B}}K(\mathcal{A},\mathcal{B}),$$

$$\Vert {\mathbf{K}}_{\nabla N\nabla M}{\mathbf{x}}_{\nabla M}-{\mathbf{y}}_{\nabla N}{\Vert }^2,$$

$$f(\mathcal{A})={\mathbf{K}}_{\mathcal{A}\nabla M}{\mathbf{x}}_{\nabla M},$$

$${\nabla }_{\mathcal{A}}f(\mathcal{A})={\mathbf{K}}_{\nabla \mathcal{A}\nabla M}{\mathbf{x}}_{\nabla M}.$$

$$\|\stackengine{\Sstackgap}{\mathbf{L}}{$\hskip 2.0pt\hat{}$}{O}{l}{F}{T}{S}\mathbf{K}\mathbf{x}-\mathbf{y}\|^{2}_{\mathbf{\Lambda}^{-1}}+\|\mathbf{x}\|^{2}_{\mathbf{K}},$$

$$\mathbf{x}={\big{[}\mathbf{K}_{MM}+(\stackengine{\Sstackgap}{\mathbf{L}}{$\hskip 2.0pt\hat{}$}{O}{l}{F}{T}{S}\mathbf{K}_{NM})^{\top}\mathbf{\Lambda}^{-1}\stackengine{\Sstackgap}{\mathbf{L}}{$\hskip 2.0pt\hat{}$}{O}{l}{F}{T}{S}\mathbf{K}_{NM}\big{]}}^{-1}(\stackengine{\Sstackgap}{\mathbf{L}}{$\hskip 2.0pt\hat{}$}{O}{l}{F}{T}{S}\mathbf{K}_{NM})^{\top}\mathbf{\Lambda}^{-1}\mathbf{y}\,.$$

$$\mathbf{y}=[\begin{array}{ccccccc}E& f_{Ax}& f_{Ay}& f_{Az}& f_{Bx}& f_{By}& f_{Bz}\end{array}{]}^{\top }.$$

$$\begin{array}{rl}{\mathbf{f}}_A& =\frac{\partial {\epsilon }_A}{\partial {\mathbf{r}}_A}+\frac{\partial {\epsilon }_B}{\partial {\mathbf{r}}_A},\\ {\mathbf{f}}_B& =\frac{\partial {\epsilon }_A}{\partial {\mathbf{r}}_B}+\frac{\partial {\epsilon }_B}{\partial {\mathbf{r}}_B}.\end{array}$$

$$\mathbf{K}=\left[\begin{array}{cc}K(\mathcal{A},\mathcal{A})& K(\mathcal{A},\mathcal{B})\\ K(\mathcal{B},\mathcal{A})& K(\mathcal{B},\mathcal{B})\end{array}\right],$$

$$\stackengine{\Sstackgap}{\mathbf{L}}{$\hskip 2.0pt\hat{}$}{O}{l}{F}{T}{S}=\left[\begin{matrix}1&1\\ \overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{A}}&\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{A}}\\ \overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{B}}&\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{B}}\\ \end{matrix}\right],$$

$$\stackengine{\Sstackgap}{\mathbf{L}}{$\hskip 2.0pt\hat{}$}{O}{l}{F}{T}{S}\mathbf{K}=\left[\begin{matrix}K({\mathcal{A}},{\mathcal{A}})+K({\mathcal{A}},{\mathcal{B}})&K({\mathcal{B}},{\mathcal{A}})+K({\mathcal{B}},{\mathcal{B}})\\ \overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{A}}K({\mathcal{A}},{\mathcal{A}})+\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{A}}K({\mathcal{B}},{\mathcal{A}})&\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{A}}K({\mathcal{A}},{\mathcal{B}})+\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{A}}K({\mathcal{B}},{\mathcal{B}})\\ \overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{B}}K({\mathcal{A}},{\mathcal{A}})+\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{B}}K({\mathcal{B}},{\mathcal{A}})&\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{B}}K({\mathcal{A}},{\mathcal{B}})+\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{B}}K({\mathcal{B}},{\mathcal{B}})\\ \end{matrix}\right]$$

$${\mathbf{K}}_{\nabla \mathcal{A}\nabla \mathcal{B}}=\left[\begin{array}{cc}{\overleftarrow{\nabla }}_{{\mathbf{r}}_A}{\overrightarrow{\nabla }}_{{\mathbf{r}}_A}K(\mathcal{A},\mathcal{A})& {\overleftarrow{\nabla }}_{{\mathbf{r}}_A}{\overrightarrow{\nabla }}_{{\mathbf{r}}_B}K(\mathcal{A},\mathcal{B})\\ {\overleftarrow{\nabla }}_{{\mathbf{r}}_B}{\overrightarrow{\nabla }}_{{\mathbf{r}}_A}K(\mathcal{B},\mathcal{A})& {\overleftarrow{\nabla }}_{{\mathbf{r}}_B}{\overrightarrow{\nabla }}_{{\mathbf{r}}_B}K(\mathcal{B},\mathcal{B})\end{array}\right].$$

$$E_{\mathrm{tot}}=\sum _i{E}^{(1)}+\sum _{i,j}{E}^{(2)}(r_{ij})+\sum _{i,j,k}{E}^{(3)}(r_{ij},r_{ik},r_{jk})+\dots$$

Full text

I Introduction

There has been a surge of activity during the last couple of years in applying machine learning methods to materials and molecular modelling problems, that was largely fuelled by the evident success of these techniques in what can loosely be called artificial intelligence. These successes have followed from the collective experience that the scientific community has gained in fitting high volumes of data with very complex functional forms that involve a large number of free parameters, while still keeping control of the regularity and thus avoiding catastrophic overfitting. In the context of molecular modelling, empirical fitting of potential energy surfaces has of course been used for many decades. Indeed it is generally held that this is the only practical way to simulate very large systems (many thousands of atoms) over long time scales (millions of time steps)Finnis (2004).

Traditionally, when fitting empirical models of atomic interactions, regularity was ensured by writing functional forms that are expressed in terms of one dimensional functions, e.g. pair potentials, spherically symmetric atomic electron densities, bond orders (as a function of number of neighbours), etc. Such functions are easy to inspect visually to ensure that they are physically and chemically meaningful, e.g. that pair potentials go to zero at large distances and are strongly repulsive at close approach, that atomic electron densities are decreasing with distance, that electron density embedding functions are convex, etc. Moreover, these natural properties are easy to build into the one dimensional functional forms or enforced as constraints in the parameter optimisation. It is widely held that employing such “physically meaningful” functional forms are key to achieving good transferability of the empirical modelsBrenner (2000).

It is also recognised, however, that the limited functional forms that can be built from these one dimensional functions ultimately limit the accuracy that these empirical models can achieve. In trying to replace them by high dimensional fits using much more flexible functional forms, two things immediately have to change. The first is the target data. When fitting only a few parameters, it is natural to demand that important observables that are deemed to be central to the scientific questions being addressed are reproduced correctly, and it is easiest to do this if they are part of the fit: e.g. melting points and other phase boundaries, radial distribution functions, etc. But in the case of very many parameters, their optimisation also takes a significant number of evaluations, and it becomes impractical to use complex observables as targets. Moreover, there is a drive towards using a “first principles” approach, i.e. that the potentials should actually reproduce the real Born-Oppenheimer potential energy surface with sufficient accuracy and therefore the scientifically relevant observables also. The hope is that this will result in transferability in the sense that a wide array of macroscopic observables will be correctly predicted without any of them being part of the fit explicitly, and therefore, the corresponding microscopic mechanisms that are also dependant on the potential energy surface will also be correct. So it is natural to take values of the potential energy, computed by the electronic structure method of choice, as the target data. The large number of free parameters can then easily be counterbalanced by a large amount of calculated target data.

The second thing that has to change is how the smooth physically meaningful behavior of the potential is controlled. It is not practical to inspect manually high dimensional functions to ensure that their predictions are physically and chemically meaningful for all possible configurations. Therefore it becomes even more important to build into the functional forms as much prior information as possible about limiting behaviour and regularity (the technical word for the kind of smoothness we are interested in). Reviewing recent work, this paper sets out an example framework for how to do this. The key goals are to create functional forms that preserve the (i) invariance of the properties over permutation of like atoms, (ii) invariance of scalar and covariance of tensorial properties with three dimensional rotations, (iii) continuity and regularity with respect to changes in atomic coordinates, including compact support of atomic interactions by including finite cutoffs.

Evidence is accumulating that strictly enforcing these physically-motivated properties is enormously beneficial, and many of the most successful machine-learning schemes for atomic-scale systems are built around symmetry arguments. One possible approach is to describe the system in terms of internal coordinates – that automatically satisfy rotational invariance – and then symmetrize explicitly the vector of representations or the functional relation between the representations and the properties. Permutationally-invariant polynomials are an example that have been very effective to model the potential energy surfaces of small molecules (see e.g. the work of Bowman and BraamsBraams and Bowman (2009)). Sorting the elements of the representation vector according to interatomic distances has also been used as a way of obtaining permutation invariance at the cost of introducing derivative discontinuities Rupp et al. (2012); Faber et al. (2015); Zhang et al. (2018). Another possibility, which we will focus on in this paper, starts from a representation of each structure in terms of atomic densities – that are naturally invariant to atom permutations – and then builds a representation that is further invariant to translations and rotations also.

Either way, once an appropriate description of each structure has been obtained, further regularization can be achieved at the level of the regression scheme. To this end, two prominent techniques are the use of artificial neural networks and kernel ridge regressionBishop (2016). We use the latter formalism here, and many further details about these techniques can be found in the rest of this volume. The kernel approach starts with the definition of a kernel function, which will be combined with a set of representative atomic configurations to construct the basis functions for the fit. It is a scalar function—at least when learning scalar quantities—with two input arguments, in the present case two atomic structures. Its value should quantify the similarity of the atomic configurations represented by its two arguments, and it can (but does not have to) be defined starting from their associated representations. The value should be largest when its two arguments are equal (or equivalent up to symmetry operations), and smallest for maximally different configurations. The degree to which the kernel is able to capture the variation of the function when varying the atomic configuration will determine how efficient the fit is. The better the correspondence, the fewer basis functions that are needed to achieve a given accuracy of fit.

II Kernel fitting

We start by giving a concise account of the kernel regression fitting approach, for more details see Refs. Bishop, 2016; Rasmussen and Williams, 2006; Schölkopf and Smola, 2002. A function defined on an atomic structure is represented as a linear sum over kernel basis functions,

[TABLE]

where the sum runs over a representative set of configurations $M$, selected from the total set $N$ of input configurations. The set of coefficients, combined into a vector $\mathbf{x}$, are determined by solving the linear system that is obtained when the available data (e.g. values of the target function evaluated for a set of structures) are substituted into Eq. (1). In the simplest case, there is one input data value corresponding to each atomic configuration. Let $\mathbf{y}$ be the vector of all available input data, and $\mathbf{K}$ be the kernel matrix with rows and columns corresponding to atomic structures, so that the element of $\mathbf{K}$ with row and column corresponding to structures ${\mathcal{A}}$ and ${\mathcal{B}}$, respectively, is $K({\mathcal{A}},{\mathcal{B}})$. The fit is then obtained by solving a linear system in the least squares sense, i.e. minimizing the quadratic loss function,

[TABLE]

The text book case is when the set of all configurations for which we have target data available is used in its entirety as the representative set (i.e. $N=M$), $\mathbf{K}$ is square, and as long as it is invertible, the optimal solution is

[TABLE]

In practice, for large data sets, using all the configurations in the data set as representatives is unnecessary. In this case, $M\subset N$, the solution is given by the pseudoinverse,

[TABLE]

where we used subscripts to emphasize the set that the vector elements correspond to, e.g. $\mathbf{y}\equiv\mathbf{y}_{N}$ is the data vector with one element for each input data structure and $\mathbf{x}\equiv\mathbf{x}_{M}$ is the vector of coefficients, one for each representative configuration. The subscripts on the kernel matrix denote array slices, i.e. $\mathbf{K}_{MN}=\mathbf{K}_{NM}^{\top}$ is the rectangular matrix whose elements correspond to the kernel values between the representative configurations and the input configurations.

Using a representative set much smaller than the total number of structures has significant advantages in terms of computational cost, often with no reduction in fitting accuracy. The training cost is dominated by computing the pseudoinverse, which scales as $O(NM^{2})$, which is linear in the size of the training data, $N$, and evaluating the model scales as $O(M)$, now independent of the size of the training data. These cost scalings are analogous to those of artificial neural networks with a fixed number of nodes.

While the above solutions are formally correct, it is widely recognised that they lead to numerical instability and overfitting, i.e. they are solutions that attempt to maximise the fit to the input data, even when this might not be desirable, which is almost always the case. At first sight, this might sound surprising, since electronic structure calculations can be made deterministic, with precise convergence behaviour in terms of its parameters, such as k-point sampling, SCF tolerance, etc. However, practical calculations are never converged to machine precision, and the resulting inconsistencies between the potential energy values for different configurations is not something that is desirable to propagate to a fitted potential energy surface. The magnitude of such inconsistencies can be easily assessed before the fit is made. Previous experienceSzlachta et al. (2014); Dragoni et al. (2018) suggests that for large databases for materials applications using plane wave density functional theory, the error due to k-point sampling are dominant, and difficult to reduce below a few meV/atom due to the associated computational cost.

In case we are fitting a potential energy surface with a representation that does not characterise the atomic positions of the whole system completely due to e.g. a finite cutoff, or some other choices made to gain computational efficiency, the fit is not expected to be exact, irrespective of the amount of input data. Sometimes, such model error can also be assessed a priori, e.g. in the case of a finite cutoff by measuring the contribution made to forces on an atom by other atoms beyond the cutoff.Bernstein et al. (2009); Deringer and Csányi (2017); Fujikake et al. (2018)

These two considerations suggest that allowing some “looseness” in the linear system might be beneficial, because it can be exploited to allow smaller linear coefficients, making the fit more regular and thus better at extrapolation. We collect the errors we expect in the fit of each target data value on the diagonal of an $N\times N$ matrix, $\mathbf{\Lambda}$. The common procedure to regularising the problem is due to TikhonovTikhonov et al. (1995). Specifically, in “kernel ridge regression” (and the equivalent “Gaussian process regression”, a Bayesian view of the same) the Tikhonov matrix is chosen to be the kernel matrix between the $M$ representative points, $\mathbf{K}_{MM}$. With highly regular (“smooth”) kernel functions, this regularisation leads to smooth fits, and the sizes of the elements of $\mathbf{\Lambda}$ control the trade-off between the accuracy of the fit and smoothness. The corresponding solutions are

[TABLE]

for the square problem, and

[TABLE]

for the rectangular problem, where we again emphasized the index sets. This solution is equivalent to minimizing

[TABLE]

which shows that the inverse of the tolerances in $\mathbf{\Lambda}$ are equivalent to regression weights on the different data points.

With the solution of the linear system in hand, the value of the fitted function for a new structure ${\mathcal{C}}$ can be written as

[TABLE]

Note that the $\mathbf{K}_{{\mathcal{C}}M}$ slice is just a vector, with elements given by the kernel between the new structure ${\mathcal{C}}$ and the structures in the representative set $M$.

II.1 Selection of a representative set

Next we describe some ways to choose the set of representative environments over which the sum in Eq. (1) is taken. This can be done by simple random sampling, but we find it advantageous to use this freedom to optimise interpolation accuracy. Among the many strategies that have been proposedHartigan and Wong (1979); Prabhakaran et al. (2012), we discuss two that have been used successfully in the context of potential energy fitting. One approach to this is to maximise the dissimilarity between the elements of the representative set. A greedy algorithm to select the configurations for the representative set is “farthest point sampling”, in which we start with a randomly selected structure, and then iteratively pick as the next structure to include the one which is farthest away from any of the structures already in the set Gonzalez (1985); Ceriotti et al. (2013); Bartók et al. (2017). The distance between two structures is measured using the “kernel metric”Schölkopf and Smola (2002), defined as

[TABLE]

This algorithm performed well for selecting molecules in regression tasks, enabling the significant reduction of the data set sizes needed to achieve a given level of accuracyDe et al. (2016).

Another technique that has been successfully used is based on matrix factorisation, which is particularly appealing when the kernel function is linear or a low order polynomial of the representation vector. Consider the matrix of feature vectors, $\mathbf{D}$, in which each row is the representation vector of an input atomic configuration, such that a linear kernel is $\mathbf{K}=\mathbf{D}\mathbf{D}^{\top}$. We are looking to select rows, many fewer than the total number, which span as much of the space as all rows span. This is a problem of matrix representation, specifically the representative set should serve as a low rank approximation of $\mathbf{K}$ and/or $\mathbf{D}$. One solution to this is called CUR matrix decompositionMahoney and Drineas (2009), which can be applied to either $\mathbf{K}$ or $\mathbf{D}$, the latter being much cheaper if the length of the representation vectors is less than the number of data points.

To determine the optimal set of representative configurations, we start with a singular value decomposition of $\mathbf{D}$,

[TABLE]

For each data point, a leverage score is calculated, essentially the weight that the top singular vectors have on that configuration.

[TABLE]

where $u_{\mathcal{A}}^{\xi}$ is the element of the $\xi$-th left singular vector that corresponds to structure ${\mathcal{A}}$. The sum runs over the first $k$ singular vectors, e.g. $k=20$ is typical. The configuration ${\mathcal{A}}$ is included in the representative set with a probability that is proportional to its leverage score, $\pi_{\mathcal{A}}$. A deterministic variant is to select one structure ${\mathcal{A}}$ at a time – the one with the highest leverage score – delete the associated row from the representation matrix and orthogonalize the remaining rows of $\mathbf{D}$ relative to it. The next data point can then be selected repeating the same procedure on the now smaller feature matrix.Imbalzano et al. (2018)

Note that in the Gaussian process literature, using a subset of the data to construct the basis is called sparsificationQuinonero-Candela and Rasmussen (2005); Snelson and Ghahramani (2005), even though the approximation relies on a low rank matrix reconstruction rather than the kernel matrix being sparse.

II.2 Linear combination of kernels

When fitting interatomic potentials for materials, a model is constructed for the atomic energy, sometimes called the “site energy”. This is both for computational efficiency and to reduce the complexity of the functional relation between structures and properties: each atomic energy is only a function of a limited number of degrees of freedom corresponding to the coordinates of the neighbouring atoms, and can therefore be evaluated independently from any other atomic energy. In fact this is the defining characteristic of an interatomic potential, in contrast to a quantum mechanical model that explicitly includes delocalised electrons. Going from atomic energies to the total energy is trivial, the latter being the sum of the former. However, going in the other direction is not unambiguous. The total energy can be calculated from a quantum mechanical model, but the atomic energies are not defined uniquely, and it becomes part of the fitting task to find the best possible decomposition of the total energy into atomic energies. Treating these two transformations on the same footing helps. Suppose we want to predict the sum of function values for two (or more) configurations. For the simple case of the sum of two energies for structures ${\mathcal{A}}$ and ${\mathcal{B}}$, the prediction is, trivially, just the sum of the individual function value predictions, e.g.

[TABLE]

If we define a new “sum-kernel” to be the sum of kernel values between a number of new configurations and the representative set, the expression for the above total energy prediction takes the same form as the prediction of the individual function values. For some set $I$ of new configurations, let

[TABLE]

where ${}^{\Sigma}\mathbf{K}_{M}$ is the vector of sum-kernel values, each element of which is the sum of the kernel between all the configurations in $I$ and a given configuration in the representative set $M$. The predicted total energy of the configurations in $I$ is then

[TABLE]

This same sum-kernel can be used to fit the model to sum data, rather than to individual function values. This is necessary in order to fit interatomic potentials for materials systems, since only total energies, and not the atomic energies themselves, are available from electronic structure calculations. At the same time, in order to enforce a finite short range in the interatomic potential, we must express the potential as an atomic energy. Using the sum-kernel, this is straightforward, the original functional form in Eq. (1) can be retained, and the we now minimise (omitting the regularisation term for clarity)

[TABLE]

where ${}^{\Sigma}\mathbf{K}$ is a matrix containing the sum-kernel values for all configurations in the input database and the representative set, and the vector $\mathbf{E}_{\mathrm{tot}}$ is the collection of corresponding total energy data.

II.3 Derivatives

The explicit analytic functional form of Eq. (1) leads to analytic derivates with respect to the atomic coordinates, e.g. forces in the case of fitting an energy. Considering for the moment the simpler case in which we are computing the derivatives of an atom-centered quantity $f({\mathcal{A}})$, we define $\nabla_{\mathcal{A}}$ as the vector of derivatives with respect to all the atomic coordinates in structure ${\mathcal{A}}$. We use the notation $\overleftarrow{\nabla{\mskip-3.0mu}}$ to indicate a derivative operator that applies to the first argument of the kernel, and $\overrightarrow{\nabla{\mskip-3.0mu}}$ to indicate a derivative that applies to the second argument. The derivatives of $f({\mathcal{A}})$ are non-zero only for atoms that belong to the structure ${\mathcal{A}}$, and are then given by differentiating Eq. (1)

[TABLE]

where we introduced the notation $\mathbf{K}_{\nabla{\mskip-3.0mu}{\mathcal{A}}M}$ to indicate the matrix that contains the derivatives of the kernels relative to all the relevant atomic coordinates. Similarly to the case of sums above, the gradient of the kernel function can also be used for fitting the model not to target values, but to gradient dataSolak et al. (2003). This is especially useful when the target represents a potential energy surface. When using typical electronic structure methods, the cost of computing the gradient with respect to all atomic positions is only a little bit more than the cost of computing the energy, but yields much more information, $3n$ pieces of data for an $n$-atom structure. There are two approaches one can take to incorporate gradient information. In the first one, used in Ref. Bartók et al., 2010 and subsequent work of that groupBartók et al. (2013); Szlachta et al. (2014); Bartók and Csányi (2015); Deringer and Csányi (2017); John and Csányi (2017); Fujikake et al. (2018); Dragoni et al. (2018); Deringer et al. (2018); Caro et al. (2018); Rowe et al. (2018); Nguyen et al. (2018), the functional form for the energy is again retained to be the same as in Eq. (1). The corresponding loss function (again without regularisation) is

[TABLE]

where $\mathbf{y}_{\nabla{\mskip-3.0mu}N}$ refers to the concatenated vector of gradients on all atoms in the set of input structures and $\mathbf{K}_{\nabla{\mskip-3.0mu}NM}$ to the corresponding matrix of kernel derivatives. The form of the solution for the coefficients is unchanged from Eq. (5) or Eq. (6) with $\mathbf{K}_{\nabla{\mskip-3.0mu}NM}$ taking the role of $\mathbf{K}_{NM}$.

In the second approach, used recently in Ref. Chmiela et al., 2017, derivatives of the kernel are the basis functions in the functional form of the fit,

[TABLE]

where $\mathbf{x}_{\nabla{\mskip-3.0mu}{\mathcal{B}}}$ contains one weight for each of the derivatives relative to the atoms in structure ${\mathcal{B}}$. The number of basis functions and corresponding coefficients is now much larger, $3nM$, for $n$-atom structures. Since the model is fitted to the derivatives, given by gradients of Eq. (18), the loss is

[TABLE]

the target properties can be computed as

[TABLE]

and their derivatives as

[TABLE]

The original motivation for this approach is apparent from Eq. (19) in which the matrix can be understood as a kernel directly between atomic forces (and in case of $M=N$, between the input data forces).

Both approaches constitute valid ways of learning a function from data representing its gradients, differing only in the choice of the kernel basis. The kernel-derivative basis functions could also be used in conjunction with a reduced representative set, and it is not yet clear which approach is better, or indeed a combination: one could choose different basis functions (kernels or their derivatives) depending on the amount and kind of data available and on the size and choice of the representative set.

II.4 Learning from linear functionals

We can combine the sum-kernel and the derivative kernel naturally, and write a single least squares problem for the coefficients in Eq. (1) that is solved to fit an interatomic potential to all available total energy, force, and virial stress data (the only condition being that the input data has to be expressible using a linear operator on function values). We define $\mathbf{y}$ as the vector with $L$ components containing all the input data: all total energies, forces and virial stress components in the training database, and $\mathbf{y}^{\prime}$ as the vector with $N$ components containing the unknown atomic energies of the $N$ atomic environments in the database, and \stackengine{\Sstackgap}{\mathbf{L}}{\hskip 2.0pt\hat{}}{O}{l}{F}{T}{S} as the linear differential operator of size $L\times N$ which connects $\mathbf{y}$ with $\mathbf{y}^{\prime}$ such that {\stackengine{\Sstackgap}{\mathbf{L}}{\hskip 2.0pt\hat{}}{O}{l}{F}{T}{S}}\mathbf{y}^{\prime}=\mathbf{y} (note that the definition of \stackengine{\Sstackgap}{\mathbf{L}}{\hskip 2.0pt\hat{}}{O}{l}{F}{T}{S} we use here is the transpose of that in Ref. Bartók and Csányi, 2015). The regularised least squares problem is now to minimise

[TABLE]

and the expression for the coefficients is given by

[TABLE]

It is instructive to write down the above matrices for the simple case when the system consists of just two atoms, $A$ and $B$, with position vectors $\mathbf{r}_{A}$, $\mathbf{r}_{B}$, target total energy $E$, and target forces $\mathbf{f}_{A}\equiv(f_{Ax},f_{Ay},f_{Az})$ and $\mathbf{f}_{B}\equiv(f_{Bx},f_{By},f_{Bz})$. The data vector is then given by

[TABLE]

The aim of the fit is to determine two unknown atomic energy functions $\varepsilon_{A}$ and $\varepsilon_{B}$ as a function of the atomic environments centered around the two atoms, ${\mathcal{A}}$ and ${\mathcal{B}}$ respectively. The total energy is their sum, $E=\varepsilon_{A}+\varepsilon_{B}$, and the forces need to include the cross terms,

[TABLE]

The representative set in this case consists of the same two atoms, so $N=M$, and the kernel matrix is square,

[TABLE]

and the linear operator \stackengine{\Sstackgap}{\mathbf{L}}{\hskip 2.0pt\hat{}}{O}{l}{F}{T}{S} is a $7\times 2$ matrix and is given by

[TABLE]

so the \stackengine{\Sstackgap}{\mathbf{L}}{\hskip 2.0pt\hat{}}{O}{l}{F}{T}{S}\mathbf{K} matrix to be substituted into equation Eq. (23) is

[TABLE]

Note that terms such as $\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{A}}K({\mathcal{B}},{\mathcal{B}})$ or $\overleftarrow{\nabla{\mskip-3.0mu}}_{\mathbf{r}_{A}}K({\mathcal{A}},{\mathcal{B}})$ are not zero because atom $A$ is present in the environment ${\mathcal{B}}$ of atom $B$, and so $K({\mathcal{B}},{\mathcal{A}})$, and also $K({\mathcal{B}},{\mathcal{B}})$, depend on $\mathbf{r}_{A}$ explicitly.

Using the approach of Ref. Chmiela et al., 2017 for the dimer, the kernel matrix is $6\times 6$, and is given by

[TABLE]

In practice it is always worth using all available data, even though once the fit is converged in the limit of infinite amount of data, the information from derivatives (forces) is the same as from energies. With finite amount of data, however, choosing the weights corresponding to energies and forces via the diagonal regulariser allows control of the fit, in the sense of its relative accuracy in reproducing energies and forces.

II.5 Learning multiple models simultaneously

Being able to fit to sums of function values has an interesting consequence. It enables in a very natural way the fitting of a model that is explicitly and a priori written as a sum of terms, each using a different kernel function, perhaps even using a different representation. That this is a good idea for potential energy functions is shown by the relative success of empirical force fields both for materials and molecules, in which the total energy is written as a sum of body-ordered terms, i.e. an atomic term, a pair potential and a three-body (angle-dependent) term, etc

[TABLE]

It is notable that while pair potentials and three body potentials using various simple parametrisations are widely used in the materials modelling literature, there are few models that take advantage of the full three dimensional flexibility of the three body term. Four-body terms are almost always restricted to one-dimensional parametrisations such as a dihedral angle. The reason for this is presumably because there is little intuition about what kinds of functional forms would be appropriate—kernel fitting avoids this problem. Such a framework was introducedBartók and Csányi (2015) and is beginning to be used for low body order model fittingDeringer and Csányi (2017); Glielmo et al. (2018); Zeni et al. (2018); Fujikake et al. (2018). Furthermore, by bringing everything together under the kernel formalism, the above expansion can also be augmented with a many-body term which enables the systematic convergence to the true Born-Oppenheimer potential energy surface, but with the many-body term having a relatively small magnitude (because the low body order terms account for most of the energy already), which helps transferability and stability.

The two-body term could be represented as a linear combination of kernels whose arguments are simply the interatomic distances, the three-body term is again a linear combination of kernels whose arguments are some representation of the geometry of three atoms, e.g. the one above using the three distances, but two distances and an angle are equally viable. The fit is then made to the target data of total energies and forces of atomic configurations, in complete analogy with Eq. (12), and now the value of the sum-kernel is the sum of pair- and triplet-kernel values between all pairs and triplets present in the two atomic configurations. A stringent test of this scheme is that in case of a target potential energy surface that is explicitly the sum of two- and three-body terms, the fit recovers these terms explicitly from just the total energies and forces Bartók and Csányi (2015).

III Density-based representations and kernels

Having summarized the algorithms that can be used to perform kernel ridge regression using atomic-scale properties and their derivatives as inputs, we now proceed to describe a framework for defining physics-based representations of local atomic environments and the kernels built from them. In kernel ridge regression, the representations do not necessarily need to be expressed explicitly, but can also be defined implicitly by means of the kernel function $K(\mathcal{A},\mathcal{B})$, that corresponds to the scalar product of representation vectors that span a (possibly infinite-dimensional) Hilbert space Schölkopf and Smola (2002). Vectors $\ket{\mathcal{A}}$ in this “reproducing kernel Hilbert space” do correspond to atomic structures, and one can write formally $K(\mathcal{A},\mathcal{B})\equiv\bra{\mathcal{A}}\ket{\mathcal{B}}$ even if the kernel might be computed without ever determining the vectors explicitly.

The reader trained in quantum mechanics will recognize an isomorphism between representations and the state vectors on one hand, and kernels and expectation values on the other. This analogy suggests that it may be beneficial to formulate atomic-scale representations using a formalism that mimics Dirac notation. Whereas in a quantum mechanical setting the physical symmetries of the problem are built into the Hamiltonian, in a machine-learning setting they are more conveniently included in the representation itself, that should be made invariant to basic symmetries such as atom labelling, rigid translations and rotations. In this section we show how starting from these intuitions one can build a very abstract description of a molecular structure that is naturally invariant with respect to the physical symmetries, based on a representation of the atom density.

Translational and rotational symmetries can be included by decomposing the structure into a collection of local environments, and by explicit symmetrization over the $SO(3)$ group. In fact, it has been recently shown Willatt et al. (2018a) how this construction leads naturally to the SOAP representation and kernel Bartók et al. (2013), and to several other popular choices of density-based representations – from Behler-Parrinello symmetry functions Behler and Parrinello (2007), to voxel density representations Kajita et al. (2017) to the binning of the pair correlation function Schütt et al. (2014) – that can be regarded as different projections of the same smooth atomic amplitude. A peculiarity of the SOAP framework is that it is formulated very naturally in terms of a kernel, that corresponds to the symmetrized overlap of atomic densities, and that it allows one to explicitly compute the representations whose scalar product constitutes the kernel function, which allows one to go back and forth between a kernel and a representation language. The atomic environmental representations can then be modified to generate non-linear kernels, as well as combined into global structural kernels. We will briefly discuss different possible approaches to the latter, either by simple linear combination of the local representations, or by a more sophisticated procedure that takes into account the most effective matching between pairs of environments in the two structures that are being compared.

III.1 A Dirac notation for structural representations

Let us introduce an abstract notation to describe atomistic structures in terms of the positions and chemical nature of the atoms that compose them Willatt et al. (2018a). Taking inspiration from Dirac notation for quantum mechanical states, we associate a ket $\ket{\mathcal{A}}$ with each configuration. Let us start with a simple example to see how such a formalism can be introduced and used. Much like in the case of quantum states, we can define a concrete representation of the ket associated with a structure in terms of positions and chemical species, e.g.

[TABLE]

where the position of each atom is represented by a smooth density $g_{i}$ (that in principle could depend on the nuclear charge and the position of atom $i$) and the kets $\ket{\alpha_{i}}$ contain the information on the nuclear charge of each atom.

The Dirac notation lends itself naturally to the definition of overlap kernels between structures, $\bra{\mathcal{A}}\ket{\mathcal{B}}$. To compute such an integral, one can use the position representation and assume that the kets associated with different elements are orthonormal:

[TABLE]

This density-based representation would not be in itself very useful, as the kernel is not invariant to relative rotations of the structures, and not even to the absolute position of the two structures in space, or their periodic representation. Nevertheless, it can be taken as the starting point to introduce many of the most successful feature representations that have been used in recent years for machine-learning of materials and molecules.

To see how, one can take inspiration from linear-scaling electronic structure methods, and the nearsightedness principle for electronic matter Yang (1991); Galli and Parrinello (1992); Goedecker (1999); Prodan and Kohn (2005). We then shift the attention from the description of complete structures to that of spherical atomic environments, that one can conveniently center on top of each atom. An atom-centered representation arises naturally from the symmetrization over the translation group of tensor products of the representation Eq. (31) Willatt et al. (2018a), and is also consistent with the atom-centered potentials that have been discussed in the previous Section as an obvious application of this framework.

We will use the notation $\ket{{\mathcal{X}_{j}}}$ to indicate an environment centered around the $j$-th atom in a structure, and express it in the position representation as

[TABLE]

where $f_{c}(r_{ij})$ is a cutoff function that restricts the environment to a spherical region centered on the atom, for the sake of computational efficiency and/or localization of the density information. The atom-centered smoothing functions are typically taken to be uniform-width Gaussians, but it would be easy to generalize the expression to include a dependency on the atomic species and/or the distance of an atom from the center of the environment, which could be used to e.g. reduce the resolution of the representation at the periphery of the environment, or adapt the smoothing length scale to each atomic species.

Note that one could also combine the density contributions from atoms of the same species into a species-dependent atomic amplitude,

[TABLE]

and then write

[TABLE]

This notation is very useful to reveal how different representations can be seen as alternative representations of the same abstract ket. For instance, one can expand the atom density in orthogonal radial functions $R_{n}(r)$ and spherical harmonics. The coefficients in such an expansion can be written as

[TABLE]

As another example, Behler-Parrinello atom-centered symmetry functions, that have been used in the construction of artificial neural network based interatomic potentials for materials Behler and Parrinello (2007); Eshet et al. (2012); Morawietz et al. (2016); Cheng et al. (2016) and molecules Smith et al. (2017) can be written by setting the basis functions to be delta distributions $g_{ij}(\mathbf{r}-\mathbf{r}_{ij})=\delta(\mathbf{r}-\mathbf{r}_{ij})$, and averaging the atom density with an appropriate pair weighting function $G_{2}$, e.g.

[TABLE]

The basis functions of the Spectral Neighbour Analysis PotentialThompson et al. (2015) also start with the same density, and expands it in hyperspherical harmonics as introduced in Ref. Bartók et al., 2010.

III.2 Smooth Overlap of Atomic Positions

It is clear that a density-based representation such as Eq. (33) is invariant to translations of the entire structure, but not to rotations that would change the orientation of the atomic neighbour amplitude. This reflects the fact that scalar products of the form $\bra{{\mathcal{X}_{j}}}\ket{{\mathcal{X}_{k}}}$ depend on the relative orientation of the environments being compared. In the Smooth Overlap of Atomic Positions (SOAP) framework, we define a symmetrized version of the overlap kernel, using the Haar integralHaar (1933) of the rotation group,

[TABLE]

where the integral is performed over all possible rotation matrices. If the base kernel is raised to the $\nu$-th power, the average preserves information on the correlations between atoms up to the $(\nu+1)$-th order Glielmo et al. (2018). As we will show below, a crucial feature of the SOAP framework is that an explicit expression for the symmetrized representation vectors $\ket{{\mathcal{X}_{j}}^{(\nu)}}$ can be given. In fact, an alternative derivation of the SOAP framework can be achieved by symmetrizing directly tensor products of the translationally invariant ket Eq. (36) Willatt et al. (2018a).

The complexity of the SOAP features is quite manageable for $\nu=1,2$, but becomes increasingly cumbersome for higher $\nu$. An effective description of higher-order interactions, that does not increase too much the complexity of the analytical evaluation of Eq. (38), can be obtained by manipulating the $\nu=2$ kernel, e.g. by taking a non-linear function of it. In practice it has been found that raising it to a power $\zeta$, and normalizing it to one

[TABLE]

is sufficient to include many-body contributions in the final kernel.

Using the Dirac notation, it is easy to see how one can give an explicit representation of the $SO(3)$ symmetrized ket for the case with $\nu=1,2$. Using a spherical harmonics expansion of $\ket{{\mathcal{X}_{j}}}$ it is very natural to perform the rotational average analytically by introducing the Wigner matrix associated with the rotation, $\bra{lm}\hat{R}\ket{l^{\prime}m^{\prime}}=\delta_{ll^{\prime}}D^{l}_{mm^{\prime}}(\hat{R})$

[TABLE]

which simplifies greatly due to the properties of the Wigner matrices. Only the term with $l=0$ survives, which makes it possible to write explicitly the $\nu=1$ symmetrized SOAP representations in terms of the spherical harmonics coefficients

[TABLE]

which corresponds to the simple kernel

[TABLE]

A position representation of the $\nu=1$ representation $\bra{r}\ket{{\mathcal{X}_{k}}^{(1)}}$ yields naturally the rotational average of $\bra{\mathbf{r}}\ket{{\mathcal{X}_{k}}}$. This can be seen by expressing $K^{(1)}({\mathcal{X}_{j}},{\mathcal{X}_{k}})$ in a position basis

[TABLE]

where we have defined the rotationally-averaged atom density

[TABLE]

which is thus closely related to the pair correlation function around the tagged atom. Similar representations have been used for machine-learning of molecules and materials Schütt et al. (2014); Faber et al. (2015), revealing once more the intimate relationships between different atom-density based representations.

The $\nu=1$ representation integrates away all angular correlations and therefore does not provide a unique representation of an environment. The representations with $\nu=2$ provide information on 3-body correlations, and can also be obtained relatively easily in closed form. The Haar integral now contains the product of two Wigner matrices. Exploiting their orthogonality relations, one obtains

[TABLE]

where the $\nu=2$ symmetrized SOAP representations read

[TABLE]

This notation corresponds to the power-spectrum components introduced in Refs. Bartók et al. (2013); De et al. (2016), $\bra{\alpha n\alpha^{\prime}n^{\prime}l}\ket{{\mathcal{X}_{j}}^{(2)}}\equiv p^{\alpha\alpha^{\prime}}_{nn^{\prime}l}({\mathcal{X}_{j}})$. Note also that, while the representation of the symmetrized kets in terms of the $nlm$ expansion is very convenient, it is not the only possibility. Similar to Eq. (44), an explicit position representation can be obtained for $\bra{\alpha\mathbf{r}_{1}\alpha^{\prime}\mathbf{r}_{2}}\ket{{\mathcal{X}_{k}}^{(2)}}$, that provides a complete representation of the 3-body rotationally-invariant correlations. The 3-body symmetry functions of the Behler-Parrinello kind can be seen as projections of this representation, similarly to the case of 2-body functions in Eq. (37).

The case of $\nu=3$ leads to an explicit representation of the ket that is proportional to the bispectrum of the environment Bartók et al. (2013)

[TABLE]

While the dimensionality of this representation makes it impractical unless somehow sparsified, it does give direct access to higher-order correlations. An interesting detail is that $\ket{{\mathcal{X}_{j}}^{(3)}}$, contrary to the $\nu=1,2$ cases, is not invariant to mirror symmetry, which makes it capable of distinguishing enantiomers.

Finally, one should note that the normalization of the kernel Eq. (39) can be achieved by normalizing the SOAP vector, so that an explicit representation of the normalized feature vector is possible. While in principle one could write out an explicit representation that yields the kernel for $\zeta>1$, it would contain an exponentially increasing number of terms. As in the case of $\ket{{\mathcal{X}_{j}}^{(3)}}$, this only makes sense if combined with a sparsification procedure.

III.3 Kernel operators and feature optimization

Provided one takes a long-range environmental cutoff, and chooses a kernel that can represent high orders of many-body interactions, a density-based representation of atomic structures should provide a complete description of any atomic structure and – given a sufficiently complete training set – predict any atomistic property with arbitrary accuracy. In practice, obviously, the accuracy of a model depends on the details of the representation, which is why different representations or kernels provide different levels of accuracy for the same training and test set Faber et al. (2017). The performance of a set of representations can be improved by modifying them so that they represent more efficiently the relations between structure and properties.

This kind of optimizations are best understood in terms of changes to the translationally-invariant environmental ket Eq. (33), and can be described, in an abstract and basis-set independent manner as a Hermitian operator acting on the ket,

[TABLE]

The most general form of this operator that makes it rotationally-invariant – so that it commutes with the rotation matrix in the definition of the SOAP kernel Eq. (38) – is readily expressed in the $\{\ket{\alpha nlm}\}$ basis Willatt et al. (2018a):

[TABLE]

While this is the most general form of the operator that is consistent with $SO(3)$ symmetry, one can use simpler forms to represent feature space transformations that can be easily understood. For instance, taking

[TABLE]

corresponds to a scaling of the smooth atom density according to the distance from the center. This kind of scaling has been shown to improve significantly the performance of SOAP Willatt et al. (2018b), as well as other density-based representations Huang and von Lilienfeld (2016); Faber et al. (2018).

Another simple form of the transformation matrix involves only the “chemical” channels

[TABLE]

This operator amounts at a change of representation for the elemental space. It is easy to see that $\bra{\alpha}\hat{U}^{\dagger}\hat{U}\ket{\alpha^{\prime}}=\kappa_{\alpha\alpha^{\prime}}$ is nothing but the “alchemical similarity matrix” between elements that has been shown to improve the accuracy of SOAP in the presence of multiple atomic species De et al. (2016); Bartók et al. (2017). What is more, by writing a low-rank approximation of $\hat{U}\approx\sum_{J\alpha}u_{J}\alpha\ket{J}\bra{\alpha}$ one can express atomic density in terms of a small number of “chemical archetypes”, improving dramatically both the accuracy and the computational cost of machine-learning models that involve more than a handful of elements Willatt et al. (2018b). Note that this transformation can be applied at the level of the translationally-invariant representation, where one can write

[TABLE]

that makes it evident how the action of this particular $\hat{U}$ operator amounts at using linear combination of atomic densities in which each specie is given weights that can be optimized by cross-validation.

III.4 $\lambda$-SOAP: Symmetry-Adapted Gaussian Process Regression

When building a machine-learning model for a tensorial property $\mathbf{T}$, one should consider that the target is not invariant under the action of a symmetry operation (e.g. a rotation) but transforms covariantly. The most effective strategy to encode the appropriate covariance properties in the model involves the decomposition of the tensor into its irreducible spherical components, i.e. combinations of the elements of the tensor that transform as the spherical harmonics of order $\lambda$ Varshalovich et al. (1988). For these irreducible components,

[TABLE]

As shown in Ref. Glielmo et al. (2017) for the case of vectors and in Ref. Grisafi et al. (2018) for tensors of arbitrary order, one has to consider a matrix-valued kernel that describes the geometric relationship between the different components of $\mathbf{T}_{\lambda}$, which can be obtained by including an additional Wigner matrix $D^{\lambda}_{\mu\mu^{\prime}}(\hat{R})$ in the Haar integral

[TABLE]

For the case with $\nu=2$ the symmetrized kets can be written explicitly based on a $\alpha nlm$ expansion of the atom density

[TABLE]

We write Eq. (55) in this form because it is somewhat symmetric, but the properties of the CG coefficients require that $m^{\prime}=m+\mu$ so the expression can be evaluated with a single sum. Furthermore, the expression evaluates to zero whenever $\left|l-l^{\prime}\right|<\lambda$, which reduces the number of elements that must be evaluated and stored, and makes it clear that Eq. (55) reduces to the scalar SOAP Eq. (46) when $\lambda=0$.

When using a linear model, each of the the symmetry-adapted representations Eq. (55) can be used to represent tensorial components that transform as $Y^{\lambda}_{\mu}$. Linearity, in this case, is necessary for preserving the symmetry properties of the $\lambda$-SOAP alpha-arxiv. A non-linear model, however, can be obtained by scaling each $\bra{\alpha nl\alpha^{\prime}n^{\prime}l^{\prime}}\ket{\mathcal{X}_{j,\lambda\mu}^{(2)}}$ by a (in principle different) non-linear function of some $\lambda=0$ representations. In the kernel language, a high-order version of the $\lambda$-SOAP kernel can be introduced with an expression analogous to Eq. (39):

[TABLE]

where $\left\|\cdot\right\|_{F}$ indicates the Frobenius norm and $\bra{{\mathcal{X}_{j}}^{(2)}}\ket{{\mathcal{X}_{k}}^{(2)}}_{\zeta-1}$ is a (scalar) SOAP kernel. This second term makes the overall kernel non-linear, without affecting the symmetry properties of the overall tensorial kernel.

III.5 Computing SOAP representations efficiently

A practical calculation of both scalar and tensorial $\nu=2$ SOAP representations $\bra{\alpha nl\alpha^{\prime}n^{\prime}l^{\prime}}\ket{\mathcal{X}_{j,\lambda\mu}^{(2)}}$ requires the evaluation of the expansion coefficients $\bra{\alpha nlm}\ket{{\mathcal{X}_{j}}}$. Let us start with the atom density written in the position representation, according to Eq. (36), and consider the case in which $\psi^{\alpha}_{\mathcal{X}}(\mathbf{r})$ is written as a superposition of spherical Gaussian functions of width $\sigma$ placed at the positions of the atoms of type $\alpha$. Then, the spherical harmonics projection in Eq. (36) can be carried out analytically, leading to:

[TABLE]

where the sum runs over all neighbouring atoms of type $\alpha$ and $\iota_{l}$ indicates a modified spherical Bessel function of the first kind. It is convenient to choose a form for the orthogonal radial basis functions $R_{n}(r)$ that makes it possible to perform the radial integration analytically.

One possible choice starts by using Gaussian type orbitals as non-orthogonal primitive functions $\tilde{R}_{k}(r)$

[TABLE]

where $\mathcal{N}_{k}$ is a normalization factor, such that $\int_{0}^{\infty}drr^{2}\tilde{R}^{2}_{k}(r)=1$. The set of Gaussian widths $\{\sigma_{k}\}$ can be chosen to span effectively the radial interval involved in the environment definition. Assuming that the smooth cutoff function approaches one at a distance $r_{\text{cut}}-\delta r_{\text{cut}}$, one could take $\sigma_{k}=(r_{\text{cut}}-\delta r_{\text{cut}})\max(\sqrt{k},1)/n_{\text{max}}$, that gives functions that are peaked at equally-spaced positions in the range between 0 and $r_{\text{cut}}-\delta r_{\text{cut}}$.

While the $\tilde{R}_{k}(r)$ are not themselves orthogonal, they can be used to write orthogonal basis functions $R_{n}(r)=\sum_{k}S^{-1/2}_{nk}\tilde{R}_{k}(r)$, where the overlap matrix $S_{kk^{\prime}}=\int\mathrm{d}rr^{2}\tilde{R}_{k}(r)\tilde{R}_{k^{\prime}}(r)$ can be computed analytically. The full decomposition of the translationally-invariant environmental ket can then be obtained without recourse to numerical integration.

Once the spherical decomposition of the atomic density has been obtained, the coefficients can be combined to give the SOAP representations of order 1 and 2. Particularly in the presence of many different chemical species, the number of components can become enormous. Ignoring for simplicity a few symmetries, and the fact that if all species do not appear in every environment it is possible to store a sparse representation nof the representation, the power spectrum contains a number of components of the order of $n_{\text{species}}^{2}n_{\text{max}}^{2}l_{\text{max}}$, which can easily reach into the tens of thousands. In the case of the tensorial $\lambda$-SOAP the number increases further to $\lambda^{2}n_{\text{species}}^{2}n_{\text{max}}^{2}l_{\text{max}}$. It is however not necessary to compute and store all of these representations: each of them, or any linear combination, is a spherical invariant (covariant) description of the environment and can be used separately as a representation. This can be exploited to reduce dramatically the computational cost and the memory footprint of a SOAP calculation, determining a low-rank approximation of the representation. One can use dimensionality reduction techniques similar to those discussed in Section II.1 to identify the most suitable reference structures. As shown in Ref. 23, both CUR decomposition and a greedy selection strategy based on farthest point sampling make it possible to reduce by more than 95% the number of SOAP representations that are needed to predict the energy of small organic molecules with chemical accuracy.

III.6 Back to the structures

Whenever one is interested in computing properties that are associated to individual atoms (for instance their NMR chemical shieldings, or the forces) one can use directly the representations corresponding to each environment, or the kernel between two environments, as the basis for a linear or non-linear regression model. As discussed in Section II, it is often the case that one is interested in using as structure labels some properties that are instead associated with the entirety of a structure – e.g. its cohesive energy, its dielectric constant, etc. In these cases a ridge regression model should be used that is based on “global” kernels between the structures, $K(\mathcal{A},\mathcal{B})$, rather than those between individual atom-centered environments. This is reflected in how the kernels between environments should be combined to give a kernel that is suitable to represent the relation between local environments and the overall property of a structure. When the target property can be seen as an additive combination of local, atom-centered contributions, the most natural (and straightforward) choice, that is consistent with Eq. (13), is

[TABLE]

It is worth stressing that in the case where the environment kernel is a linear kernel based on SOAP representations, this sum-kernel can be written in terms of a global representation associated with the entire structure,

[TABLE]

where we introduced

[TABLE]

An alternative way to combine the information from individual environments in a symmetrized global kernel corresponds to averaging the Fourier coefficients of each environment,

[TABLE]

and then taking the Haar integral of the resulting sum. For instance, for $\nu=2$,

[TABLE]

The form Eq. (59) is more general, and one can readily introduce non-linear kernels such as $\bra{{\mathcal{X}_{j}}^{(\nu)}}\ket{{\mathcal{X}_{k}}^{(\nu)}}_{\zeta}$ for which an explicit expression for the representations would be too cumbersome. Eq. (59) also suggests that the combination of environment kernels could be generalized by introducing a weighting matrix

[TABLE]

One could for instance determine the importance of each environment within a structure, and set $W_{jk}(\mathcal{A},\mathcal{B})=w_{j}(\mathcal{A})w_{k}(\mathcal{B})$. Alternatively, one can use techniques from optimal transport theory Cuturi (2013) to define an entropy-regularized matching (REMatch) procedure De et al. (2016), in which $W_{jk}$ is a doubly stochastic matrix that matches the most similar environments in the two structures, disregarding the environmental kernels between very dissimilar environments

[TABLE]

where $d^{2}$ indicates the kernel-induced squared distance Eq. (9). The parameter $\gamma$ weights the entropy regularization and makes it possible to interpolate between strict matching of the most similar pairs of environments ($\gamma\rightarrow 0$) to an average kernel that weights all pairs equally ($\gamma\rightarrow\infty$). Although this construction complicates considerably the combination of local kernels, it provides a strategy to introduce an element of non-locality in the comparison between structures. Given the cost of computing the REMatch kernel, and the fact that it prevents using some sparsification strategies that act at the level of individual environments, this method should be used when the target property is expected to exhibit very strong non-additive behavior, e.g. when just one portion of the system is involved – for instance when determining the activity of a drug molecule, a problem for which REMatch has been shown to improve dramatically the accuracy of the ML model Bartók et al. (2017).

III.7 Multi-kernel learning

We have shown that SOAP representations can be seen as just one possible embodiement of a general class of rotationally-symmetrized density-based representations, that also encompasses other popular representations for atomic-scale machine learning, and that can be tuned to a great extent, e.g. by changing the way different components are weighted. The fact that different representations can be computed within the same formalism does not imply they are fully equivalent: each expression or kernel emphasizes different components of the structure/property relations. For instance, kernels with varying radial scaling or cutoff distance focus the machine-learning model on short, mid or long-range interactions. It is then natural to consider whether a better overall model can be constructed by combining representations that are associated with different cutoff distances, or different levels of body order expansions. This can be achieved by a weighted combination of kernels of the form

[TABLE]

where each $K_{\aleph}$ corresponds to a distinct model.

This is equivalent to an additive model for a property, similar to the construction of an atom-centered decomposition of the total energy in Eq. (12). In this case, instead, the property $y$ associated with each structure is written as the sum of contributions $y_{\aleph}({\mathcal{A}})$ that are associated with the various kernels $K_{\aleph}$

[TABLE]

where $x_{{\mathcal{B}}}$ are the kernel regression weights for each of the representative structures ${\mathcal{B}}$. The weights $w_{\aleph}$ correspond to the estimated contribution that each model will give to the final property, and can be obtained by cross-validation, or by physical intuition. For instance, in the case of multiple radial cutoffs, it is found that much smaller weights should be associated with long-range kernels, consistent with the fact that distant interactions contribute a small (although often physically relevant) contribution to the total energy Bartók et al. (2017). It should also be noted that, provided that the representations corresponding to the kernels are linearly independent, Eq. (66) effectively corresponds to a feature space of increased dimensionality, obtained by concatenating the representations that are – implicitly or explicitly – associated with each kernel.

Conclusions

We have laid out a mathematical framework, based on the concept of the atomic density, for building representations of atomic environments that preserve the geometric symmetries, and chemically sensible limits. Coupled with kernel regression, this allows the fitting of complex models of physical properties on the atomic scale, both scalars like interatomic potentials (force fields), and tensors such as multipole moments and quantum mechanical operators. We discuss in general terms how kernel regression can be extended to include a sparse selection of reference structures, and to predict and learn from linear functionals of the target property. To leverage the many formal similarities between kernel regression and quantum mechanics, we use a Dirac bra-ket notation to formulate the main results concerning the SOAP representations. This notation also helps making apparent the relationship between SOAP representations and other popular density-based approaches to represent atomic structures. The framework can be extended and tuned in many different ways to incorporate insight about the relations between properties, structures and representations. With physical principles such as symmetry and nearsightedness of interactions at its core, we believe this formulation is ideally suited to provide a unified framework to machine learn atomic-scale properties.