Machine Learning Unifies the Modelling of Materials and Molecules

TL;DR

This paper introduces a machine learning framework that unifies the modeling of materials and molecules, accurately predicting atomic properties and stability across different chemical systems using a local environment description and Bayesian learning.

Contribution

It presents a novel, universal machine learning approach that captures quantum effects and predicts stability and activity with high accuracy, advancing atomistic modeling.

Findings

Predicts silicon surface reconstructions with quantum accuracy

Accurately assesses molecular stability across classes

Distinguishes active/inactive protein ligands with >99% reliability

Abstract

Determining the stability of molecules and condensed phases is the cornerstone of atomistic modelling, underpinning our understanding of chemical and materials properties and transformations. Here we show that a machine learning model, based on a local description of chemical environments and Bayesian statistical learning, provides a unified framework to predict atomic-scale properties. It captures the quantum mechanical effects governing the complex surface reconstructions of silicon, predicts the stability of different classes of molecules with chemical accuracy, and distinguishes active and inactive protein ligands with more than 99% reliability. The universality and the systematic nature of our framework provides new insight into the potential energy surface of materials and molecules.

Figures

Equations

$$\bar{y}(x)=\sum _i{w}_{i}K(x,x_i),$$

$$p(y^{\ast }|\mathbf{y})\propto p(y^{\ast }\&\mathbf{y})=\mathcal{N}(\bar{y},{\sigma }^{\ast 2})$$

$$K(A,B)=\sum _{i\in A,j\in B}{P}_{ij}k\left(x_i,x_j\right).$$

Full text

Machine Learning Unifies the

Modelling

of Materials and Molecules

Albert P. Bartók,1 Sandip De,2,3 Carl Poelking,4

Noam Bernstein,5 James Kermode,6 Gábor Csányi,7 Michele Ceriotti3,∗

1Scientific Computing Department, Science and Technology Facilities Council,

Rutherford Appleton Laboratory, Oxfordshire OX11 0QX, United Kingdom

2National Center for Computational Design and Discovery

of Novel Materials (MARVEL)

3Laboratory of Computational Science and Modelling,

Institute of Materials, EPFL, Lausanne, Switzerland

4Department of Chemistry,

University of Cambridge, Cambridge CB2 1EW, United Kingdom

5Center for Materials Physics and Technology,

U.S. Naval Research Laboratory, Washington, DC 20375, USA

6Warwick Centre for Predictive Modelling,

School of Engineering, University of Warwick, Coventry CV4 7AL, United Kingdom

7Engineering Laboratory, University of Cambridge

∗To whom correspondence should be addressed; E-mail: michele.ceriotti@epfl.ch

Determining the stability of molecules and condensed phases is the cornerstone of atomistic modelling, underpinning our understanding of chemical and materials properties and transformations. Here we show that a machine-learning model, based on a local description of chemical environments and Bayesian statistical learning, provides a unified framework to predict atomic-scale properties. It captures the quantum mechanical effects governing the complex surface reconstructions of silicon, predicts the stability of different classes of molecules with chemical accuracy, and distinguishes active and inactive protein ligands with more than 99% reliability. The universality and the systematic nature of our framework provides new insight into the potential energy surface of materials and molecules.

RESEARCH SUMMARY: Statistical learning based on a local representation of atomic structures provides a universal model of chemical stability

1 Introduction

Calculating the energies of molecules and of condensed-phase structures is fundamental to predicting the behavior of matter at the atomic scale, and a formidable challenge. Reliably assessing the relative stability of different compounds, and of different phases of the same material, requires the evaluation of the energy of a given three-dimensional assembly of atoms with an accuracy comparable with the thermal energy ($\sim$0.5 kcal/mol at room temperature), which is a small fraction of the energy of a chemical bond (up to $\sim$230 kcal/mol for the \ceN2 molecule).

Quantum mechanics is a universal framework that can deliver this level of accuracy. By solving the Schrödinger equation, the electronic structure of materials and molecules can in principle be computed, and from it all ground-state properties and excitations follow. The prohibitive computational cost of exact solutions at the level of electronic-structure theory lead to the development of many approximate techniques that address different classes of systems. Coupled-cluster theory (CC) (?) for molecules, and density-functional theory (DFT) (?, ?, ?), for the condensed phase have been particularly successful and can typically deliver the levels of accuracy required to address a plethora of important scientific questions. The computational cost of these electronic structure models is nevertheless still significant, limiting their routine application in practice to dozens of atoms in the case of CC and hundreds in the case of DFT.

To go further, explicit electronic structure calculation has to be avoided, and we have to predict the energy corresponding to an atomic configuration directly. While such empirical potential methods (force fields) are indeed much less expensive, their predictions to date have been qualitative at best. Moreover, the number of distinct approaches have rapidly multiplied – in the struggle for accuracy at low cost, generality is invariably sacrificed. Recently, machine-learning approaches have started to be applied to designing interatomic potentials that interpolate electronic-structure data as opposed to using parametric functional forms tuned to match experimental or calculated observables. While there have been several hints that this approach can achieve the accuracy of DFT at a fraction of the cost (?, ?, ?, ?, ?, ?, ?), little effort has been put into recovering the generality of quantum mechanics: atomic and molecular descriptors, and learning strategies have been optimized for different classes of problems, and in particular efforts for materials and for chemistry have been rather disconnected. Here we show that the combination of Gaussian process regression (?) with a local descriptor of atomic neighbor environments that is general and systematic can re-unite the modelling of hard matter and molecules, consistently achieving predictive accuracy.

This lays the foundations for a universal reactive force field that can recover the accuracy of the Schrödinger equation at negligible cost and – thanks to the locality of the model – leads to an intuitive understanding of the stability and the interactions between molecules. By showing that we can accurately classify active and inactive protein ligands we provide evidence that this framework can be extended to capture more complex, non-local properties as well.

2 Results

2.1 The reconstructions of silicon surfaces

Due to its early technological relevance to the semiconductor industry and simple bulk structure, Si has traditionally been one of the archetypical tests for new computational approaches to materials modelling  (?, ?, ?, ?, ?, ?). In spite of the fact that its bulk properties can be captured reasonably well by simple empirical potentials, its surfaces display remarkably complex reconstructions, whose stability is governed by a subtle balance of elastic properties and quantum mechanical effects, such as the Jahn-Teller distortion that determines a tilt of dimers on Si(100). The determination of the dimer-adatom-stacking fault (DAS) $7\times 7$ reconstruction of Si(111) as the most stable among several similar structures was the culmination of a concerted effort of experiment and modelling including early scanning tunnelling microscopy  (?), and was also a triumph for DFT  (?).

As shown in Figure 1, empirical potentials incorrectly predict the un-reconstructed $1\times 1$ to be a lower energy configuration, and fail to predict the $7\times 7$ as the lowest energy structure even from among the DAS reconstructions. Up to now, the only models that could capture these effects included electronic structure information, at least on the tight binding level (or its approximation as a bond order potential). We trained a SOAP-GAP model on a database of configurations from short ab initio molecular dynamics trajectories of small unit cells (including the $3\times 3$ reconstruction, but not those with larger unit cells; for details, see SI). This model correctly describes a broad array of standard bulk and defected material properties within a wide range of pressures and temperatures, as well as properties that depend on transition state energetics such as the generalized stacking fault surfaces shown in the SI. A striking illustration of the power of this model is the quantitative description of both the tilt of the (100) dimers and the ordering of the (111) reconstructions – without explicitly considering the quantum mechanical electron density.

Nevertheless, even this model is based on a training dataset which is a result of ad hoc (if well informed) choices. The Bayesian GPR framework tells us how to improve the model. The predicted error $\sigma^{*}$, shown as the color scale in Fig. 1b, can be used to identify new configurations likely to be usefully added to the training set. The adatoms on the surface have the highest error, and once we included small surface unit cells with adatoms, the ML model came much closer to its target.

2.2 Coupled-cluster energies for 130k Molecules

Molecular properties exhibit distinctly different challenges than bulk materials, from the combinatorial number of stable configurations, to the presence of collective quantum mechanical and electrostatic phenomena such as aromaticity, charge transfer and hydrogen bonding. At the same time, many relevant scientific questions involve predicting the energetics of stable conformers, which is a less-complex problem than obtaining a reactive potential. Following early indication of success on a small dataset  (?, ?), here we start our investigation using the GDB9 dataset that contains about 134,000 small organic molecules whose geometries have been optimized at the level of DFT, and that has been used in many of the pioneering studies of machine learning for molecules  (?, ?). Accurate models have been reported, however, only when predicting DFT-energies using as inputs geometries that have already been optimized at the DFT level – which makes the exercise insightful  (?) but does not constitute an alternative to doing the DFT calculation.

Figure 2a demonstrates that the GPR framework using the very same SOAP descriptors can be used to obtain useful predictions of the chemical energy of a molecule (the atomization energy) on this heterogeneous chemical dataset. DFT methods give very good equilibrium geometries, and are often used as a stepping stone to evaluate energies at the “gold standard” level of CC theory (CCSD(T)). They have also been shown to constitute an excellent baseline reference towards higher levels of theory  (?). Indeed, a SOAP-GAP model can use DFT inputs and only 500 training points to predict CCSD(T) atomization energies with an error below the symbolic threshold of 1 kcal/mol. The error drops to less than 0.2 kcal/mol when training on 15% of the GDB9.

DFT calculations for the largest molecules in GDB9 can nowadays be performed in a few hours, which is still impractical if one wanted to perform high-throughput molecular screening on millions of candidates. Instead, we can use the inexpensive semi-empirical PM7 model (taking around a second to compute a typical GDB9 molecule) to obtain an approximate relaxed geometry, and build a model to bridge the gap between geometries and energies (?). With a training set of  20,000 structures, the model predicts CCSD(T) energies with 1 kcal/mol accuracy using only the PM7 geometry and energy as input.

In order to achieve this level of accuracy it is however crucial to use this information judiciously. The quality of PM7 training points, as quantified by the root-mean square difference (RMSD) $d$ between PM7 and DFT geometries, varies significantly across the GDB9. In keeping with the Bayesian spirit of the ML framework, we set the diagonal variance $\propto\exp(d^{2}/\lambda^{2})$ corresponding to the prior information that structures with a larger RMSD between the two methods may be affected by a larger uncertainty. Even though we do not use RMSD information on the test set, the effect of down-weighting information from the training structures for which PM7 gives inaccurate geometries is to reduce the prediction error by more than 40%.

The strategy used to select training structures also has a significant impact on the reliability of the model. Fig. 2b shows a sketch-map  (?) of the structure of the GDB9 dataset based on the kernel-induced metric, demonstrating the inhomogeneity of the density of configurations. Random selection of reference structures leaves large portions of the space unrepresented, which results in a very heavy-tailed distribution of errors (see SI). We find that selecting the training set sequentially using a greedy algorithm that picks the next farthest data point to be included (farthest point sampling, FPS) gives more uniform sampling of the database, dramatically reducing the fraction of large errors, especially in the peripheral regions of the dataset (Fig. 2c and d), leading to a more resilient ML model. It should be noted that this comes at the price of a small degradation of the performance as measured by the commonly used mean absolute error (MAE), due to the fact that densely populated regions do not get any preferential sampling.

In order to test the “extrapolative power”, or transferability of the SOAP-GAP framework we then applied the GDB9-trained model for $\Delta_{\text{DFT-CC}}$ to the prediction of the energetics of larger molecules, and considered $\sim 850$ conformers of the dipeptides obtained from two natural amino acids, aspartic acid and glutamic acid  (?). Although GDB9 does not explicitly contain information on the relative energies of conformers of the same molecule, we could predict the CCSD(T) corrections to the DFT atomization energies with an error of 0.45 kcal/mol, a 100-fold reduction compared to the intrinsic error of DFT.

It is worth stressing that, within the scope of the SOAP-GAP framework, there is considerable room for improvement of the accuracy. Using the same SOAP parameters that we adopted for the GDB9 model for the benchmark task of learning DFT energies using DFT geometries as inputs, we could obtain a mean absolute error of 0.40 kcal/mol in the smaller QM7b dataset  (?). As discussed in the SI, using an “alchemical kernel”  (?) to include correlations between different species allowed us to further reduce that error to 0.33 kcal/mol. A “multi-scale” kernel (a sum of SOAP kernels each with a different radial cutoff parameter) that combines information from different length scales allows one to reach an accuracy of 0.26 kcal/mol (or alternatively, to reach 1 kcal/mol accuracy with fewer than 1000 FPS training points) – both results being considerably superior to existing methods that have been demonstrated on similar datasets. The same multi-scale kernel also improves significantly the performance for GDB9, allowing us to reach 1 kcal/mol with just 5000 reference energies, and as little as 0.18 kcal/mol with 75000 structures. Given that SOAP-GAP allows naturally to both predict and learn from derivatives of the potential (i.e. forces), the doors are open for building models that can describe local fluctuations and/or chemical reactivity by extending the training set to non-equilibrium configurations – as we demonstrated already for the silicon force field here, and previously for other elemental materials.

2.3 The stability of molecular conformers

To reduce even further the prediction error on new molecules, we could include a larger set of training points from the GDB9. It is clear from the learning curve in Fig. 2a that the ML model is still far from its saturation point. For the benchmark DFT learning exercise we could attain an error smaller than 0.28 kcal/mol using 100k training points, which can be improved even further using a more complex multi-scale kernel (see SI). An alternative is to train a specialized model that aims to obtain accurate predictions of the relative energies of a set of similar molecules. As an example of this approach, we considered a set of 208 conformers of glucose, whose relative stability has been recently assessed with a large set of electronic-structure methods  (?). Fig. 3a shows that as few as 20 reference configurations are sufficient to evaluate the corrections to semiempirical energies that are needed to reach 1 kcal/mol accuracy relative to complete-basis-set CCSD(T) energies, or to reach 0.2-0.4 kcal/mol error when using different flavors of DFT as a baseline.

2.4 Receptor ligand binding

The accurate prediction of molecular energies opens up the possibility of computing a vast array of more complex thermodynamic properties, using the SOAP-GAP model as the underlying energy engine in molecular dynamics simulation. However, the generality of the SOAP kernel for describing chemical environments also allows directly attacking different classes of scientific questions – e.g. sidestepping not only the evaluation of electronic structure, but also the cost of demanding free-energy calculations, making instead a direct connection to experimental observations. As a demonstration of the potential of this approach, we investigated the problem of ligand-receptor binding. We used data from the DUD-E (Directory of Useful Decoys, Enhanced)  (?), a highly curated set of ligand-receptor pairs taken from the ChEMBL database, enriched with property-matched decoys  (?). These decoys resemble the individual ligands in terms of atomic composition, molecular weight, and physicochemical properties, but are structurally distinct in that they do not bind to the protein receptor.

We trained a Kernel-Support-Vector-Machine (Kernel-SVM) (?, ?) for each of the 102 receptors listed in the DUD-E, to predict whether or not each candidate molecule binds to the corresponding polypeptide. We used an equal but varying number $n_{\textrm{train}}$ of ligands and decoys (up to 120) for each receptor, using the SOAP kernel as before to represent the similarity between atomic environments. Here however we chose the matrix $\mathbf{P}$ in eq. (3) corresponding to an optimal permutation matching (“MATCH”-SOAP) rather than a uniform average  (?). Predictions are collected over the remaining compounds and the results are averaged over different subsets used for training.

The receiver-operating characteristic (ROC), shown in Fig. 4, describes the trade-off between the rate of true positives $p(+|+)$ versus false positives $p(+|-)$ as the decision threshold of the SVM is varied. The area under the ROC curve (AUC) is a widely used performance measure of binary classifiers, in a loose sense giving the fraction of correctly classified items. A SOAP-based SVM trained on just 20 examples can predict receptor ligand binding with a typical accuracy of 95%, which goes up to 98% when 60 training examples are used, and 99% when using a FPS training set selection strategy – significantly surpassing the present state-of-the-art (?, ?, ?). The model is so reliable that its failures are highly suggestive of inconsistencies in the underlying data. The dashed line in Fig. 4a corresponds to receptor FGFR1 and shows no predictive capability. Further investigation uncovered data corruption in the DUD-E dataset, with identical ligands labelled both as active and inactive. Using an earlier version of the database  (?) shows no such anomaly, giving an AUC of 0.99 for FGFR1.

3 Discussion

Machine learning is often regarded – and criticized – as the quintessentially naïve inductive approach to science. In many cases, however, one can extract some intuition and insight from a critical look at the behavior of a machine-learning model.

Fitting the difference between levels of electronic structure theory gives an indication of how smooth and localized, and therefore easy for SOAP-GAP to learn, are the corrections that are added by increasingly expensive methods. For instance, hybrid DFT methods are considerably more demanding than plain “generalized-gradient approximation” DFT, and indeed show a considerably smaller baseline variance to high-end quantum chemistry methods. However, the error of the corresponding SOAP-GAP model is almost the same for the two classes of DFT, which indicates that exact-exchange corrections to DFT are particularly short ranged, and therefore easy to learn with local kernel methods. Thanks to the additive nature of the average-kernel SOAP kernel, it is also possible to decompose the energy difference between methods into atom-centered contributions (Fig. 3b). The discrepancy between DFT and semiempirical methods appears to involve large terms with opposite sign (positive for carbon atoms, negative for aliphatic hydrogens), that partially cancel out. Exact exchange plays an important role in determining the energetics of the ring and open chain forms  (?), and indeed the discrepancy between PBE and PBE0 is localized mostly on the aldehyde/hemiacetal group, as well as, to a lesser extent, on the H-bonded O atoms. The smaller corrections between CC methods and hybrid functionals show less evident patterns, as one would expect when the corrections involve correlation energy.

Long-range non-additive components to the energy are expected for any system with electrostatic interactions – and could be treated, for instance, by machine-learning the local charges and dielectric response terms  (?), and then feeding them into established models of electrostatics and dispersion.

However for elemental materials and the small molecules we consider here an additive energy model can be improved simply by increasing the kernel range, $r_{c}$. Looking at the dependence of the learning curves on the cutoff for the GDB9 (see SI), we can observe the trade-off between the completeness of the representation and its extrapolative power  (?). For small training set sizes, a very short cutoff of 2 Å and the averaged molecular kernel give the best performance, but then saturates at about 2 kcal/mol. Longer cutoffs give initially worse performance, because the input space is larger, but the learning rate deteriorates more slowly; at 20,000 training structures, $r_{c}=3$ Å yields the best performance. Given that the SOAP kernel gives a complete description  (?) of each environment up to $r_{c}$, we can infer from these observations the relationship between the length and energy scales of physical interactions (see SI). For a DFT model, considering interactions up to 2 Å is optimal if one is content to capture physical interactions with an energy scale of the order of 2.5 kcal/mol. When learning corrections to electron correlation, $\Delta_{\text{DFT-CC}}$, most of the short-range information is already included in the DFT baseline, and so length scales up to and above 3 Å become relevant already for $n_{\text{train}}<20,000$, allowing an accuracy of less that 0.2 kcal/mol to be reached.

In contrast, the case of ligand binding predictions poses a significant challenge to an additive energy model already at the small molecule scale. Ligand binding is typically mediated by electro-negative/positive or polarizable groups located in “strategic” locations within the ligand molecule, which additionally must satisfy a set of steric constraints in order to fit into the binding pocket of the receptor. Capturing these spatial correlations of the molecular structure is a prerequisite to accurately predicting whether or not a given molecule binds to a receptor. This is demonstrated by the unsatisfactory performance of a classifier based on an averaged combination of atomic SOAP kernels (see Fig. 4b). By combining the atomic SOAP kernels using an “environment matching” procedure, one can introduce a degree of non-locality – because now environments in the two molecules must be matched pairwise, rather than in an averaged sense. Thus, the relative performance of different kernel combination strategies give a sense of whether the global property of a molecule can result from averages over different parts of the system, or whether a very particular spatial distribution of molecular features is at play.

A striking demonstration of inferring structure-property relations from a ML model is given in Fig. 4b-c, where the SOAP classifier is used to identify binding moieties (“warheads”) for each of the receptors. To this end, we formally project the SVM decision function $z$ onto individual atoms of a test compound associating to each a “binding score” (see SI). Red and yellow regions of the isosurface plots denote moieties that are expected to promote binding. For decoys, no consistent patterns are resolved. The identified warheads are largely conserved across the set of ligands – in fact, by investigating the position of the crystal ligand inside the binding pocket of the adenosine receptor A2 (b), we can confirm that a positive binding field is indeed assigned to those molecular fragments that localize in the pocket of the receptor. Scanning through the large set of ligands in the dataset (see SI), it is also clear that the six-membered ring and its amine group, fused with the adjacent five-membered ring, are the most prominent among the actives. Finally, note that regions of the active ligands colored in blue (as in Fig. 4c) could serve as target locations for lead optimization, e.g., to improve receptor affinity and selectivity.

The consistent success of the SOAP-GAP framework across materials, molecules and biological systems shows that it is possible to sidestep the explicit electronic structure and free energy calculation and determine the direct relation between molecular geometry and stability. This already enables useful predictions to be made in many problems, and there is a lot of scope for further development – e.g. by using a deep-learning approach, by developing multi-scale kernels to treat long range interactions, using active learning strategies (?), or by fine tuning the assumed correlations between the contributions of different chemical elements, as discussed in the SI. We believe that the exceptional performance of the SOAP-GAP framework we demonstrated stems from its general, mathematically rigorous approach to the problem of representing local chemical environments. Building on this local representation allowed us to capture even more complex, non-local properties.

4 Materials and Methods

Gaussian process regression (GPR) is a Bayesian machine learning framework (?) which is also formally equivalent to another machine learning method, Kernel Ridge Regression (KRR). Both are based on a kernel function $K(x,x^{\prime})$ that acts as a similarity measure between inputs $x$ and $x^{\prime}$. Data points close in the metric space induced by the kernel are expected to correspond to the values $y$ and $y^{\prime}$ of the function one is trying to approximate. Given a set of training structures $x_{i}$ and the associated properties $y_{i}$, the prediction of the property for a new structure $x$ can be written as

[TABLE]

which is a linear fit using the kernel function as a basis, evaluated at the locations of the prior observations. The optimal setting of the weight vector is ${\bf w}=({\bf K}+\sigma^{2}{\bf I})^{-1}{\bf y}$, where $\sigma$ is the Tichonov regularization parameter. In the framework of GPR, which takes as its prior a multivariate normal distribution with the kernel as its covariance, Eq. (1) represents the mean, $\bar{y}$, of the posterior distribution

[TABLE]

which now also provides an estimate of the error of the prediction, $\sigma^{*}$. The regularization parameter $\sigma$ corresponds to the expected deviation of the observations from the underlying model due to statistical or systematic errors. Within GPR it is also easy to obtain generalizations for observations that are not of the function values, but linear functionals thereof (sums, derivatives). Low-rank (sparse) approximations of the kernel matrix are straightforward and help reduce the computational burden of the matrix inversion in computing the weight vector (?).

The efficacy of machine learning methods critically depends on developing an appropriate kernel, or equivalently, on identifying relevant features in the input space that are used to compare data items. In the context of materials modelling, the input space of all possible molecules and solids is vast. We can drastically reduce the learning task by focusing on local atomic environments instead, and using a kernel between local environments as a building block.

We use the Smooth Overlap of Atomic Positions (SOAP) kernel, which is the overlap integral of the neighbor density within a finite cutoff $r_{c}$, smoothed by a Gaussian with a length scale governed by the interatomic spacing, and finally integrated over all 3D rotations and normalized. This kernel is equivalent to the scalar product of the spherical power spectra of the neighbor density (?), which therefore constitutes a chemical descriptor of the neighbor environment. Both the kernel and the descriptor respect all physical symmetries (rotations, translations, permutations), are smooth functions of atomic coordinates and can be refined at will to provide a complete description of each environment.

To construct a kernel $K$ between two molecules (or periodic structures) A and B from the SOAP kernel $k$ we average over all possible pairs of environments,

[TABLE]

As shown in the SI, choosing $P_{ij}=\frac{1}{N_{A}N_{B}}$ for fitting the energy per atom is equivalent to defining it as a sum of atomic energy contributions (i.e. an interatomic potential), with the atomic energy function being a GPR fit using the SOAP kernel as its basis. Given that the available observations are total energies and their derivatives with respect to atoms (forces), the learning machine provides us with the optimal decomposition of the quantum mechanical total energy into atomic contributions. In keeping with the nomenclature of the recent literature, we call a GPR model of the atomistic potential energy surface a “Gaussian Approximation Potential” (GAP), and a “SOAP-GAP model” is one which uses the SOAP kernel.

Other choices of $P$ are possible and will make sense for various applications. For example, setting $P$ to be the permutation matrix that maximizes the value of $K$ corresponds to the “best match” assignment between constituent atoms in the two structures that are compared - which can be computed in polynomial time by formulating the task as an optimal assignment problem  (?). It is possible to smoothly interpolate between the average and best match kernels using an entropy-regularized Wasserstein distance  (?) construction.

5 Acknowledgments

ABP was supported by a Leverhulme Early Career Fellowship and the Isaac Newton Trust until 2016, ABP also acknowledges support from CCP-NC funded by EPSRC (EP/M022501/1). SD was supported by the NCCR MARVEL, funded by the Swiss National Science Foundation. MC acknowledges funding by the European Research Council under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 677013-HBMAP). CP and JRK were supported by the EU grant “NOMAD”, grant no. 676580. JRK acknowledges support from the EPSRC under grants EP/L014742/1 and EP/P002188/1. GC acknowledges support from EPSRC grants EP/L014742/1, EP/J010847/1 and EP/J022012/1. The work of NB was funded by the Office of Naval Research through the U. S. Naval Research Laboratory’s core basic research program. Computations were performed at the Argonne Leadership Computing Facility under contract DE-AC02-06CH11357, the High Performance Computing Service at Cambridge University, computing resources provided by STFC Scientific Computing Department’s SCARF cluster, and also ARCHER under the “UKCP” EPSRC grants EP/K013564/1 and EP/P022561/1.

Data Availability All data needed to evaluate the conclusions in the paper are present in the paper, in the cited references and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Competing Interests: The authors declare that they have no competing interests.

Author contributions: APB, SD, GC and MC performed and analysed calculations on molecular databases. CP, GC and MC performed and analysed drug binding predictions. APB, NB, JRK and GC performed and analysed calculations on silicon surfaces. All the authors contributed to the writing of the manuscript.

References and Notes

A. Szabo, N. S. Ostlund, Modern Quantum Chemistry (Dover Publications, 2012).

• 2.

R. M. Martin, Electronic structure: basic theory and practical methods (Cambridge University Press, 2004).

• 3.

P. Hohenberg, W. Kohn, Phys. Rev. 136, B864 (1964).

• 4.

W. Kohn, L. J. Sham, Phys. Rev. 140, A1133 (1965).

• 5.

J. Behler, M. Parrinello, Phys. Rev. Lett. 98, 146401 (2007).

• 6.

A. P. Bartók, M. C. Payne, R. Kondor, G. Csányi, Phys. Rev. Lett. 104, 136403 (2010).

• 7.

J. C. Snyder, M. Rupp, K. Hansen, K.-R. Müller, K. Burke, Phys. Rev. Lett. 108, 253002 (2012).

• 8.

G. Montavon, et al., New Journal of Physics 15, 095003 (2013).

• 9.

F. A. Faber, A. Lindmaa, O. A. von Lilienfeld, R. Armiento, Phys. Rev. Lett. 117, 135502 (2016).

• 10.

A. V. Shapeev, Multiscale Modeling & Simulation 14, 1153 (2016).

• 11.

J. S. Smith, O. Isayev, A. E. Roitberg, Chem. Sci. 8, 3192 (2017).

• 12.

C. E. Rasmussen, C. K. I. Williams, Gaussian processes for machine learning (MIT Press, 2006).

• 13.

R. Wolkow, Phys. Rev. Lett. 68, 2636 (1992).

• 14.

G. Binnig, H. Rohrer, C. Gerber, E. Weibel, Phys. Rev. Lett. 50, 120 (1983).

• 15.

R. Car, M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).

• 16.

P. Rinke, A. Janotti, M. Scheffler, C. G. V. de Walle, Physical Review Letters 102 (2009).

• 17.

A. J. Williamson, J. C. Grossman, R. Q. Hood, A. Puzder, G. Galli, Physical Review Letters 89 (2002).

• 18.

J. Behler, R. Martonák, D. Donadio, M. Parrinello, Phys. Rev. Lett. 100, 185501 (2008).

• 19.

K. D. Brommer, M. Needels, B. Larson, J. D. Joannopoulos, Phys. Rev. Lett. 68, 1355 (1992).

• 20.

S. De, A. P. Bartók, G. Csányi, M. Ceriotti, Phys. Chem. Chem. Phys. 18, 13754 (2016).

• 21.

R. Ramakrishnan, P. O. Dral, M. Rupp, O. A. von Lilienfeld, Scientific Data 1, 1 (2014).

• 22.

R. Ramakrishnan, P. O. Dral, M. Rupp, O. A. von Lilienfeld, J. Chem. Theory Comput. 11, 2087 (2015).

• 23.

K. T. Schütt, F. Arbabzadah, S. Chmiela, K. R. Müller, A. Tkatchenko, Nature Comm. 8, 13890 (2017).

• 24.

M. Ceriotti, G. A. Tribello, M. Parrinello, Proc. Natl. Acad. Sci. USA 108, 13023 (2011).

• 25.

M. Ropo, M. Schneider, C. Baldauf, V. Blum, Scientific Data 3, 160009 (2016).

• 26.

M. Marianski, A. Supady, T. Ingram, M. Schneider, C. Baldauf, J. Chem. Theory Comput. 12, 6157 (2016).

• 27.

http://www.libatoms.org/dude-soap/.

• 28.

V.-P. Jaakola, et al., Science 322, 1211 (2008).

• 29.

M. M. Mysinger, M. Carchia, J. J. Irwin, B. K. Shoichet, Journal of Medicinal Chemistry 55, 6582 (2012).

• 30.

N. Lagarde, J.-F. Zagury, M. Montes, Journal of Chemical Information and Modeling 55, 1297 (2015).

• 31.

B. Schölkopf, A. Smola, K.-R. Müller, Neural Computation 10, 1299 (1998).

• 32.

B. Schölkopf, A. J. Smola, Learning with kernels: support vector machines, regularization, optimization, and beyond, Adaptive computation and machine learning (MIT Press, Cambridge, Mass, 2002).

• 33.

P. Skoda, D. Hoksza (IEEE, 2016), pp. 1220–1227.

• 34.

A. A. Lee, M. P. Brenner, L. J. Colwell, Proceedings of the National Academy of Sciences 113, 13564 (2016).

• 35.

I. Wallach, M. Dzamba, A. Heifets, CoRR abs/1510.02855 (2015).

• 36.

N. Huang, B. K. Shoichet, J. J. Irwin, Journal of Medicinal Chemistry 49, 6789 (2006).

• 37.

N. Artrith, T. Morawietz, J. Behler, Phys. Rev. B 83, 153101 (2011).

• 38.

B. Huang, O. A. von Lilienfeld, J. Chem. Phys. 145 (2016).

• 39.

A. P. Bartók, R. Kondor, G. Csányi, Phys. Rev. B 87, 184115 (2013).

• 40.

Z. Li, J. R. Kermode, A. De Vita, Phys. Rev. Lett. 114, 096405 (2015).

• 41.

J. Quiñonero-Candela, C. E. Rasmussen, Journal of Machine Learning Research 6, 1939 (2005).

• 42.

H. W. Kuhn, Naval Research Logistics Quarterly 2, 83 (1955).

• 43.

M. Cuturi, Advances in Neural Information Processing Systems 26, C. J. C. Burges, L. Bottou, M. Welling, Z. Ghahramani, K. Q. Weinberger, eds. (Curran Associates, Inc., 2013), pp. 2292–2300.

• 44.

W. J. Szlachta, A. P. Bartók, G. Csányi, Physical Review B 90, 104108 (2014).

• 45.

V. L. Deringer, G. Csányi, Physical Review B 95, 094203 (2017).

• 46.

J. P. Perdew, et al., Physical Review B 46, 6671 (1992).

• 47.

S. J. Clark, et al., Z. Kristall. 220, 567 (2005).

• 48.

A. Van Duin, S. Dasgupta, F. Lorant, W. Goddard III, Journal of Physical Chemistry A 105, 9396 (2001).

• 49.

M. J. Buehler, A. C. T. van Duin, W. A. Goddard, Physical Review Letters 96, 095505 (2006).

• 50.

T. Lenosky, et al., Modelling and Simulation in Materials Science and Engineering 8, 825 (2000).

• 51.

J. Tersoff, Physical Review B 38, 9902 (1988).

• 52.

F. H. Stillinger, T. A. Weber, Physical Review B 31, 5262 (1985).

• 53.

S. D. Solares, et al., Langmuir 21, 12404 (2005).

• 54.

J. Sadowski, J. Gasteiger, G. Klebe, Journal of Chemical Information and Computer Sciences 34, 1000 (1994).

• 55.

N. M. OBoyle, et al., Journal of Cheminformatics 3, 33 (2011).

• 56.

J. M. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman, D. A. Case, Journal of Computational Chemistry 25, 1157 (2004).

• 57.

J. J. P. Stewart, Journal of Molecular Modeling 19, 1 (2012).

• 58.

J. J. P. Stewart, MOPAC 2016, http://openmopac.net (2016).

• 59.

J. Frisch, M, et al., Gaussian 09, Revision D.01 (2013). Gaussian Inc. Wallingford CT 2013.

• 60.

R. Ramakrishnan, P. O. Dral, M. Rupp, O. A. von Lilienfeld, Scientific Data 1, 1 (2014).

• 61.

A. D. Becke, Journal of Chemical Physics 98, 5648 (1993).

• 62.

P. J. Stephens, F. J. Devlin, C. F. Chabalowski, M. J. Frisch, Journal of Physical Chemistry 98, 11623 (1994).

• 63.

L. A. Curtiss, K. Raghavachari, P. C. Redfern, V. Rassolov, J. A. Pople, The Journal of Chemical Physics 109, 7764 (1998).

• 64.

H.-J. Werner, et al., Molpro, version 2012.1, a package of ab initio programs (2012). See https://www.molpro.net.

• 65.

F. A. Faber, et al., J. Chem. Theory Comput. p. acs.jctc.7b00577 (2017).

• 66.

H. Huo, M. Rupp, arXiv preprint arXiv:1704.06439 (2017).

• 67.

R. Krishnan, J. S. Binkley, R. Seeger, J. A. Pople, The Journal of Chemical Physics 72, 650 (1980).

• 68.

M. Ceriotti, G. A. Tribello, M. Parrinello, J. Chem. Theory Comput. 9, 1521 (2013).

• 69.

C. Poelking, SOAPXX, https://github.com/capoe/soapxx (2017).

See pages 1 of si See pages 2 of si See pages 3 of si See pages 4 of si See pages 5 of si See pages 6 of si See pages 7 of si See pages 8 of si See pages 9 of si