Machine learning prediction errors better than DFT accuracy
Felix A. Faber, Luke Hutchison, Bing Huang, Justin Gilmer and, Samuel S. Schoenholz, George E. Dahl, Oriol Vinyals, Steven Kearnes, and Patrick F. Riley, O. Anatole von Lilienfeld
TL;DR
This study demonstrates that machine learning models can predict molecular properties with errors smaller than those of DFT calculations, potentially surpassing hybrid DFT accuracy when high-quality data is available.
Contribution
The paper systematically compares various regressors and molecular representations, showing ML models can outperform DFT in predicting electronic properties of molecules.
Findings
ML predictions deviate less than DFT from experimental data
Prediction errors are close to chemical accuracy
ML models could surpass DFT with better data
Abstract
We investigate the impact of choosing regressors and molecular representations for the construction of fast machine learning (ML) models of thirteen electronic ground-state properties of organic molecules. The performance of each regressor/representation/property combination is assessed using learning curves which report out-of-sample errors as a function of training set size with up to 117k distinct molecules. Molecular structures and properties at hybrid density functional theory (DFT) level of theory used for training and testing come from the QM9 database [Ramakrishnan et al, {\em Scientific Data} {\bf 1} 140022 (2014)] and include dipole moment, polarizability, HOMO/LUMO energies and gap, electronic spatial extent, zero point vibrational energy, enthalpies and free energies of atomization, heat capacity and the highest fundamental vibrational frequency. Various…
| Feature | Description |
| Atom type | H, C, N, O, F (one-hot). |
| Chirality | R or S (one-hot or null). |
| Formal charge | Integer electronic charge. |
| Ring sizes | For each ring size (–), the number of rings that include this atom. |
| Hybridization | , , or (one-hot or null). |
| Hydrogen bonding | Whether this atom is a hydrogen bond donor and/or acceptor (binary values). |
| Aromaticity | Whether this atom is part of an aromatic system. |
| Feature | Description |
| Bond type | Single, double, triple, or aromatic (one-hot or null). |
| Graph distance | For each distance (–), whether the shortest path between the atoms in the pair is less than or equal to that number of bonds (binary values). |
| Same ring | Whether the atoms in the pair are in the same ring. |
| Spatial distance | The euclidean distance between the two atoms. |
| NMMAE | ||||||||||||
| arb. u. | ||||||||||||
| Mean | ||||||||||||
| MAD | ||||||||||||
| Target | ||||||||||||
| DFT | - | |||||||||||
| EN | CM | |||||||||||
| BOB | ||||||||||||
| BAML | ||||||||||||
| ECFP4 | ||||||||||||
| HDAD | ||||||||||||
| HD | ||||||||||||
| MARAD | ||||||||||||
| Mean | ||||||||||||
| BR | CM | |||||||||||
| BOB | ||||||||||||
| BAML | ||||||||||||
| ECFP4 | ||||||||||||
| HDAD | ||||||||||||
| HD | ||||||||||||
| MARAD | ||||||||||||
| Mean | ||||||||||||
| RF | CM | |||||||||||
| BOB | ||||||||||||
| BAML | ||||||||||||
| ECFP4 | ||||||||||||
| HDAD | ||||||||||||
| HD | ||||||||||||
| MARAD | ||||||||||||
| Mean | ||||||||||||
| KRR | CM | |||||||||||
| BOB | ||||||||||||
| BAML | ||||||||||||
| ECFP4 | ||||||||||||
| HDAD | ||||||||||||
| HD | ||||||||||||
| MARAD | ||||||||||||
| Mean | ||||||||||||
| GG | MG | |||||||||||
| GC | MG |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Machine learning prediction errors better than DFT accuracy
Felix A. Faber
Institute of Physical Chemistry and National Center for Computational Design and Discovery of Novel Materials, Department of Chemistry, University of Basel, Klingelbergstrasse 80, CH-4056 Basel, Switzerland
Authors contributed equally
Luke Hutchison
Google, 1600 Amphitheatre Parkway, Mountain View, CA, US - 94043 CA
Authors contributed equally
Bing Huang
Institute of Physical Chemistry and National Center for Computational Design and Discovery of Novel Materials, Department of Chemistry, University of Basel, Klingelbergstrasse 80, CH-4056 Basel, Switzerland
Justin Gilmer
Google, 1600 Amphitheatre Parkway, Mountain View, CA, US - 94043 CA
Samuel S. Schoenholz
Google, 1600 Amphitheatre Parkway, Mountain View, CA, US - 94043 CA
George E. Dahl
Google, 1600 Amphitheatre Parkway, Mountain View, CA, US - 94043 CA
Oriol Vinyals
Google, 5 New Street Square, London EC4A 3TW, UK
Steven Kearnes
Patrick F. Riley
Google, 1600 Amphitheatre Parkway, Mountain View, CA, US - 94043 CA
O. Anatole von Lilienfeld
Institute of Physical Chemistry and National Center for Computational Design and Discovery of Novel Materials, Department of Chemistry, University of Basel, Klingelbergstrasse 80, CH-4056 Basel, Switzerland
Abstract
We investigate the impact of choosing regressors and molecular representations for the construction of fast machine learning (ML) models of thirteen electronic ground-state properties of organic molecules. The performance of each regressor/representation/property combination is assessed using learning curves which report out-of-sample errors as a function of training set size with up to 117k distinct molecules. Molecular structures and properties at hybrid density functional theory (DFT) level of theory used for training and testing come from the QM9 database [Ramakrishnan et al, Scientific Data 1 140022 (2014)] and include dipole moment, polarizability, HOMO/LUMO energies and gap, electronic spatial extent, zero point vibrational energy, enthalpies and free energies of atomization, heat capacity and the highest fundamental vibrational frequency. Various representations from the literature have been studied (Coulomb matrix, bag of bonds, BAML and ECFP4, molecular graphs (MG)), as well as newly developed distribution based variants including histograms of distances (HD), and angles (HDA/MARAD), and dihedrals (HDAD). Regressors include linear models (Bayesian ridge regression (BR) and linear regression with elastic net regularization (EN)), random forest (RF), kernel ridge regression (KRR) and two types of neural networks, graph convolutions (GC) and gated graph networks (GG). Out-of sample errors are strongly dependent on the choice of representation and regressor and molecular property. Electronic properties are typically best accounted for by MG and GC, while energetic properties are better described by HDAD and KRR. The specific combinations with the lowest out-of-sample errors in the 117k training set size limit are dipole moment (MG/GC), static polarizability (MG/GG), HOMO/LUMO eigenvalue and gap (MG/GC), electronic spatial extent (BOB/KRR), zero point vibrational energy (HDAD/KRR), (free) energies and enthalpies of atomization (HDAD/KRR), heat capacity at room temperature (HDAD/KRR), and highest fundamental vibrational frequency (BAML/RF). We present numerical evidence that ML model predictions deviate from DFT less than DFT deviates from experiment for all properties. Furthermore, our out-of-sample prediction errors with respect to hybrid DFT reference are on par with, or close to, chemical accuracy. Our findings suggest that ML models could be more accurate than hybrid DFT if explicitly electron correlated quantum (or experimental) data was available.
††preprint: APS/123-QED
I Introduction
Due to its favorable trade-off between accuracy and computational cost, Density Functional Theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) is the work-horse of quantum chemistry Burke (2012)—despite its well known shortcomings regarding spin-states, van der Waals interactions, and chemical reactions Koch and Holthausen (2002); Cohen et al. (2012). Failures to predict reaction profiles are particularly worrisome Plata and Singleton (2015), and recent analysis casts even more doubts on the usefulness of DFT functionals obtained through parameter fitting Medvedev et al. (2017). The prospect of universal and computationally much more efficient machine learning (ML) models, trained on data from experiments or generated at higher levels of electronic structure theory such as post-Hartree Fock or quantum Monte Carlo (e.g. exemplified in Ref. Ramakrishnan et al. (2015)), therefore represents an appealing alternative to DFT. Not surprisingly, a lot of recent effort has been devoted to developing ever more accurate ML models of properties of molecular and condensed phase systems.
Several ML studies have already been published using a data set called QM9 Ramakrishnan et al. (2014), consisting of molecular quantum properties for the 134k smallest organic molecules containing up to 9 heavy atoms (C, O, N, or F; not counting H) in the GDB-17 universe Ruddigkeit et al. (2012). Some of these studies have developed or used representations we consider in this work, such as BAML (Bonds, angles, machine learning) Huang and von Lilienfeld (2016), bag of bonds (BoB) Hansen et al. (2015); Ramakrishnan and von Lilienfeld (2015) and the Coulomb matrix (CM) Rupp et al. (2012); Ramakrishnan and von Lilienfeld (2015). Atomic variants of the CM have also been proposed and tested on QM9 Barker et al. (2016). Other representations have also been benchmarked on QM9 (or QM7 which is a smaller but similar data set), such as Fourier series of radial distance distributions von Lilienfeld et al. (2015), motifs Huan et al. (2015), the smooth overlap of atomic positions (SOAP) Bartók et al. (2013) in combination with regularized entropy match De et al. (2016), constant size descriptors based on connectivity and encoded distance distributions Collins et al. (2016). Ramakrishnan et al. (2015) introduced a -ML approach, where the difference between properties calculated at coarse/accurate quantum level of theories is being modeled. Furthermore, neural network models, as well as deep tensor neural networks, have recently been proposed and tested on the same or similar data sets Smith et al. (2017); Schütt et al. (2017). Dral et al. (2015) use such data to machine learn optimal molecule specific parameters for the OM2 Weber and Thiel (2000) semi-empirical method, and orthogonalization tests are benchmarked in Ref. Dral et al. (2016).
However, limited work has yet been done in systematically assessing various methods and properties on large sets of the exact same chemicals Hansen et al. (2013). In order to unequivocally establish if ML has the potential to replace DFT for the screening of properties, one has to demonstrate that ML test errors are systematically lower than estimated DFT accuracies for all the properties available. This study accomplishes that through a large scale assessment of unprecedented scale: (i) In order to approximate large training set sizes , we included 13 quantum properties from up to 117k molecules in training (90% of QM9). (ii) We tested multiple regressors (Bayesian ridge regression (BR), linear regression with elastic net regularization (EN), random forest (RF), kernel ridge regression (KRR), neural network (NN) models graph convolutions (GC) Kearnes et al. (2016) and gated graphs (GG) Li et al. (2016)) and (iii) multiple representations including BAML, BOB, CM, extended connectivitiy fingerprints (ECFP4), histograms of distance, angle, and dihedral (HDAD), molecular alchemical radial angular distribution (MARAD), and molecular graphs (MG). (iv) We investigated all combinations of regressors and representations, except for MG/NN which was exclusively used together because GC and GG depend fundamentally on the input representation being a graph instead of a flat feature vector.
The best models for the various properties are: dipole moment (MG/GC), static polarizability (MG/GG), HOMO/LUMO eigenvalue and gap (MG/GC), electronic spatial extent (BOB/KRR), zero point vibrational energy (HDAD/KRR), atomization energy at 0 Kelvin (HDAD/KRR), atomization energy at room temperature (HDAD/KRR), enthalpy of atomization at room temperature (HDAD/KRR), atomization of free energy at room temperature (HDAD/KRR), heat capacity at room temperature (HDAD/KRR), and the highest fundamental vibrational frequency (BAML/RF). For training set size of 117k (90 of data set) we have found the additional out-of-sample error added by machine learning to be lower or as good as DFT errors relative to experiment for all properties, and that chemical accuracy (See table 3) is reached, or in sight.
This paper is organized as follows: First we will briefly describe the methods, including data set, model validation protocols, representations, and regressors. In section III, we present the results and discuss them, and section IV concludes the paper.
II Method
II.1 Data set
We have used the QM9 data set consisting of 134k drug-like organic molecules Ramakrishnan et al. (2014). Molecules in the data set consist of H, C, O, N and F, and contain up to 9 heavy atoms. For each molecule several properties, calculated at DFT level of theory (B3LYP/6-31G(2df,p)), were included. We used: Norm of dipole moment (Debye), norm of static polarizability (Bohr3), HOMO eigenvalue (eV), LUMO eigenvalue (eV), HOMO-LUMO gap (eV), electronic spatial extent (Bohr2), zero point vibrational energy ZPVE (eV), atomization energy at 0 Kelvin (eV), atomization energy at room temperature (eV), enthalpy of atomization at room temperature (eV), atomization of free energy at room temperature (eV), heat capacity at room temperature (cal/mol/K), and the highest fundamental vibrational frequency (cm*-1*). For energies of atomization (, , and ) all models yield very similar errors: We will only discuss for for the remainder. The 3053 molecules specified in Ref. Ramakrishnan et al. (2014) which failed SMILES consistency tests were excluded from our study, as well as two linear molecules, leaving 131k molecules.
II.2 Model validation
Starting from the 131k molecules in QM9 after removing the 3k molecules (see above) we have created a number of train-validation-test splits. We have splitted the data set into test and non-test sets and varied the percentage of data in test set to explore the effect of amount of data in error rates. Inside the non-test set, we have performed 10 fold cross validation for hyperparameter optimization. That is, for each model 90% (the training set) of the non-test set is used for training and 10% (the validation set) is used for hyperparameter selection. For each test/non-test split, we have trained 10 models on different subsets of the non-test set, and we report the mean error on the test set across those 10 models. Note that the non-test set will be referred to as training set in the results section in order to simplify discussion.
In terms of CPU investments necessary for training the respective models we note that EN/BR, RF/KRR, and GC/GG required minutes, hours, and multiple days, respectively. Using GPUs could dramatically reduce such timings.
II.3 DFT errors
To place the quality of our prediction errors in the right context, experimental accuracy estimates of DFT become desirable. Here, we summarize literature results comparing DFT to experiments for the various properties we study. Where data is available, the corresponding deviation of DFT from experiment is listed in Table 3, alongside our ML prediction errors (vide infra). More specifically, rough hybrid DFT error estimates for dipole moment and polarizability have been obtained from Refs. Koch and Holthausen (2002); Hickey and Rowley (2014).
Frontier molecular orbital energies (HOMO, LUMO and HOMO-LUMO gap) can not be measured directly. However, for the exact (yet unknown) exchange-correlation potential, HOMO and LUMO eigenvalues correspond to the negative of the vertical ionization potential (IP) and electron affinity (EA), respectively Stowasser and Hoffmann (1999). In order to quantify the accuracy of DFT (B3LYP/6-31G(2df,p)) based frontier eigenvalues, we have carried out additional corresponding calculations for a set of 28 organic molecules with up to 9 heavy atoms (not counting hydrogen) for which experimental (adiabatic) IPs and EAs were reported for NIST nis (2016), or in Ref. Parr and Yang (1989) (see also Table S2 in supplementary materials). Ignoring the (possibly severe) difference between vertical and adiabatic values, we thus approximate DFT errors of HOMO/LUMO energies by their respective deviation from experimentally measured IPs and EAs. Correspondlingly, we estimate the HOMO-LUMO gap error by comparison to the experimental IP/EA gap.
DFT RMSE estimates of ZPVE and (the highest fundamental vibrational frequency) were published in Ref. Sinha et al. (2004) for a set of 41 organic molecules, with up to 6 heavy atoms (not counting hydrogen) and calculated using B3LYP/cc-pVTZ (assuming a normal distribution, we have calculated corresponding MAEs using a scaling factor of 0.8Geary (1935)). DFT errors estimates of energies of atomization have been taken from Ref. Curtiss et al. (1997). The errors were estimated from 69 organic molecules (22 hydrocarbons and 47 substituted hydrocarbons), with up to 6 heavy atoms (not counting hydrogens), and calculated using B3LYP/6-311+G(3df,2p). Deviation of DFT (at the B3LYP/6-311g** level of theory) from experimental heat capacities were reported by DeTar (2007) who obtained errors of 16 organic molecules, with up to 8 heavy atoms (not counting hydrogens). No satisfactory experimental data for electronic spacial extent has been found in the literature.
Note, however, that one should be cautious when referring to these errors: Strictly speaking they can not be compared since different basis sets, molecules, and experiments were used. Nevertheless, we feel that the quoted errors provide meaningful guidance as to what one can expect from DFT for each property.
II.4 Representations
II.4.1 CM and BoB
The Coulomb matrix (CM) representationRupp et al. (2012) is a square atom by atom matrix, where off diagonal elements are the Coulomb nuclear repulsion terms between atom pairs. The diagonal elements approximate the electronic potential energy of the free atoms. Atom indices in the CM are sorted by the norm of each atom’s row (or column). The Bag of Bonds (BoB)Hansen et al. (2015) representation uses exclusively CM elements, grouping them for different atom pairs into different bags, and sorting them within each bag by their relative magnitude.
II.4.2 BAML
The recently introduced BAML (Bonds, angles, machine learning) representation can be viewed as a many-body extension of BoBHuang and von Lilienfeld (2016). All pairwise nuclear repulsions are replaced by Morse/Lennard-Jones potentials for bonded/non-bonded atoms respectively. Furthermore, three- and four-body interactions between covalently bonded atoms are included using angular and torsional terms, respectively. Parameters and functional forms are based on the universal force field (UFF)Rappe et al. (1992).
II.4.3 ECFP4
Extended Connectivity Fingerprints Rogers and Hahn (2010) (ECFP4) are a common representation of molecules in cheminformatics based studies. They are particularly popular for drug discovery Besnard et al. (2012); Lounkine et al. (2012); Huigens III et al. (2013). The basic idea, typical also for other cheminformatics descriptors Todeschini and Consonni (2009) (e.g. the signature descriptor Faulon et al. (2003); Visco et al. (2002)) is to represent a molecule as the set of subgraphs up to a fixed diameter (here we use ECFP4, which is a max diameter of 4 bonds). To produce a fixed length vector, the subgraphs can be hashed such that every subgraph sets one bit in the fixed length vector to 1. In this work, we use a fixed length vector of size 1024. Note that ECFP4 is based solely on the molecular graph specifying all covalent bonds, e.g. as encoded by SMILES strings.
II.4.4 MARAD
Molecular atomic radial angular distribution (MARAD) is a radial distribution function (RDF) based representation. MARAD is an extension of an atomic representation ARAD, consisting of three RDFs for distances between atom pairs, and parallel and orthogonal projections of distances in atom triplets, respectively. Distances between two molecules can be evaluated analytically. Unfortunately, most regressors evaluated in this work, such as BR, EN and RF, do not rely on inner products and distances between representations. We remedy this issue by projecting MARAD onto bins in order to work with all regressors (apart for GG and GC which use MG exclusively). Further technical details regarding MARAD are specified in the Supplementary materials sup .
II.4.5 HD, HDA, and HDAD
BOB, BAML and MARAD rely on computing functions for given interatomic distances, and/or angles, and/or torsions, and then either project that value onto discrete bins, or sort the values. As a straightforward alternative, we also investigated representations which account directly from pairwise distances, triple-wise angles, and quad-wise dihedral angles through manually generated bins in histograms. The resulting representations in increasing interatomic many-body order are called HD (Histogram of distances), HDA (Histogram of distances and angles), and HDAD (Histogram of distances, angles and dihedral angles). For any given molecule, one iterates through each atom , producing a set of distances, angle and dihedral angle features for .
Distance features were produced by measuring the distance between and (for ) for each element pair. The distance features were assigned a label incorporating the atomic symbols of and sorted alphabetically (with H last), e.g. if was a carbon atom and was a nitrogen atom, the distance feature for the atom pair would be labeled C-N. These labels will be used to group all features with the same label into a histogram and allow us to only count each pair of atoms once.
Angle features were produced by taking the principal angles formed by the two vectors spanning from each atom to every subset of 2 of its 3 nearest atoms, and . The angle features were labeled by the element type of , followed by the alphabetically sorted element types (Except for hydrogens, which were listed last) of and . The example where is a Carbon atom, a Hydrogen atom, a Nitrogen would be assigned the label C-N-H.
Dihedral angle features were produced by taking the principal angles between two planes. We take as the origin, and for each of the four nearest neighbors in turn, labeling the neighbor atom , and forming a vector . Then all subsets of the remaining three out of four nearest neighbors of are chosen, and labeled as and . This third and fourth atom respectively form two triangular faces when paired with : and . The dihedral angle between the two triangular faces was calculated. These dihedral angle features were labeled with the atomic symbol for , followed by the atomic symbols for , and , sorted alphabetically, with the exception that hydrogens were listed last, e.g. C-C-N-H.
The features from all molecules have been aggregated for each label type to generate a histograms for each label type. Fig. 1 exemplifies this for C-N distances, C-C-C angles, and C-C-C-O dihedrals for the entire QM9 data set. Certain typical molecular characteristics can be recognized upon mere inspection. For example, the CN histogram displays a strong and isolated peak between 1.1 and 1.5 Å, corresponding to occurrences of single, double, and triple bonds. For distances above 2 Å, peaks at typical radii of second and third covalent bonding shells around N can be recognized at 2.6 Å and 3.9 Å. Also C-C-C angles can be easily interpreted: The peak close to zero and Rad corresponds to geometries where three atoms are part of a linear (alkyne, or nitrile) type of motif. The broad and largest peak corresponds to 120 and 109 degrees, typically observed in and hybridized atoms.
The morphology of each histogram has then been examined to identify apparent peaks and troughs, motivated by the idea that peaks indicate structural commonalities among molecules. Bin centers have been placed at each significant local minimum and maximum (Shown as vertical lines in Fig. 1). Values at 15-25 bin centers have been chosen as a representation for each label type. All bin center values are provided in the Supplementary Material. For each molecule, the collection of features has subsequently been rendered into a fixed-size representation, producing one vector component for each bin center, within each label type. This has been accomplished follwoing a two-step process. (i) Binning and interpolation: Each feature value is projected on the two nearest bins. The relative amount projected on each bin uses linear projection between the two bins. For example: A feature with value which lies between two bins placed at and respectively, contributes and to the first and second bin respectively. (ii) Reduction: The collection of contributions within each bin of each molecule’s feature vector is condensed to a single value by summing all contributions.
II.4.6 Molecular Graphs
We have investigated several neural network models which are based on molecular graphs (MG) as representation. The inputs are real-valued vectors associated with each atom and with each pair of atoms. More specifically, we have used the featurization described in Kearnes et al. (2016) with the removal of partial charge and the addition of Euclidean distances to the pair feature vectors. All elements of the feature vector are described in Tables 1 and 2.
The featurization process was unsuccessful for a small number of molecules (367) because of conversion failures from geometry to rational SMILES string when using OpenBabel O’Boyle et al. (2011) or RDKit Landrum (2014), and were excluded from all results using the molecule graph features.
Note that within a previous draft of this study Faber et al. (2017), we reported biased results for GC/GG models due to use of Mulliken partial charges within the MG representation. All MG results presented herewithin have been obtained without any Mulliken charges in the representation. Model hyper parameters for the GC model, however, still correspond to the previously obtained hyper parameter search.
II.5 Regressors
For all methods, we first standardized the property values so that all properties have zero mean and unit standard deviation.
II.5.1 Kernel Ridge Regression
KRR Müller et al. (2001); Schölkopf and Smola (2002); Vovk (2013); Hastie et al. (2011) is a type of regression with regularization Hoerl et al. (2000) which uses kernel functions as basis set. A property of a query molecule is predicted by a sum of weighted kernel functions between and all molecules in the training set,
[TABLE]
where are regression coefficients, obtained by minimizing the Euclidean distance between the estimated and the reference property of all molecules in the training set. We used Laplacian and Gaussian kernels as implemented by scikit-learn Pedregosa et al. (2011) for all representations.
The level of noise in our data is very low so strong regularization is not necessary. We therefore set the regularization parameter to . Kernel widths were chosen by screening values on a base-2 logarithmic grid for the 10 percent training set (from 0.25 to 8192 for Gaussian kernel and 0.1 to 16384 for Laplacian kernel). In order to simplify the width screening, prior to learning all feature vectors were normalized (scaling the input vector by the mean norm across the training set) by the Euclidean norm for the Gaussian kernel and the Manhattan norm for the Laplacian kernel.
II.5.2 Bayesian Ridge Regression
We use BR Neal (1996) as is implemented in scikit-learn Pedregosa et al. (2011). BR is a linear model with a penalty on the coefficients. Unlike Ridge Regression where the strength of that penalty is a regularizing hyperparameter which must be set, in Bayesian Ridge Regression the optimal regularizer is estimated from the data.
II.5.3 Elastic Net
Also EN Zou and Hastie (2005) is a linear model. Unlike BR, the penalty on the weights is a mix of and terms. In addition to the regularizing hyperparameter for the weight penalty, Elastic net has an additional hyperparameter l1_ratio to control the relative strength of the and weight penalties. We used the scikit-learn Pedregosa et al. (2011) implementation and set . We then did a hyperparmaeter search on regularizing parameter in a base 10 logarithmic grid from to .
II.5.4 Random Forest
We use RF Breiman (2001) as implemented in scikit-learn Pedregosa et al. (2011). RF regressors produce a value by averaging many individual decision trees fitted on randomly resampled sets of the training data. Each node in each decision tree is a threshold of one input feature. Early experiments did not reveal strong differences in performance based on the number of trees used, once a minimal number was reached. We have used 120 trees for all regressions.
II.5.5 Graph Convolutions
We have used the GC model as described in Kearnes et al. (2016), with several structural modifications and optimized hyperparameters. The graph convolution model is built on the concepts of “atom” layers (one real vector associated with each atom) and “pair” layers (one real vector associated with each pair of atoms). The graph convolution architecture defines operations to transform atom and pair layers to new atom and pair layers. There are three structural changes to the model used herewithin when compared to the one described in Kearnes et al. (2016). We describe these briefly here with details in the Supplementary Material. First, we have removed the “Pair order invariance” property by simplifying the () transformation. Since the model only uses the atom layer for the molecule level features, pair order invariance is not needed. Second, we have used the Euclidean distance between atoms. In the () transformation, we divide the value from the convolution step by a series of distance exponentials. If the original convolution for an atom pair with distance produces a vector , we concatenate the vectors , , , , and to produce the transformed value for the pair . Third, we have followed other work on neural networks based on chemical graphs Duvenaud et al. (2015). Inspired by fingerprinting like Extended Connectivity Fingerprints Rogers and Hahn (2010), based on summing softmax operations to convert a real valued vector to a sparse vector and sum those sparse vectors across all the atoms, we have used the same operation here along with a simple sum across the atoms to produce molecule level features from the top atom layer. We have found that this works as well or better than the Gaussian histograms first used in GC Kearnes et al. (2016). To optimize the network, we have searched the hyperparameter space using Gaussian Process Bandit Optimization Desautels et al. (2014) as implemented by HyperTune hyp (2016). The hyperparameter search has been based on the evaluation of the validation set for a single fold of the data. Further details including parameters, and search ranges chosen for this paper are listed in the Supplementary materials sup .
II.5.6 Gated Graph Neural Networks
We have used the GG Neural Networks model (GG) as described in Li et al. (2016). Similar to the GC model, it is a deep neural network whose input is a set of node features , and an adjacency matrix with entries in a discrete set to indicate different edge types. It has internal hidden representations for each node in the graph of dimension which it updates for steps of computation. Its output is invariant to all graph isomorphisms, meaning the order of the nodes presented to the model does not matter. We have implemented the GG model as described in the original work with some additional preprocessing of the input. To include distance information we bin atomic distances into 10 bins, [0, 2], [2,2.5], [2.5,3], [3,3.5], [3.5,4], [4,4.5], [4.5,5], [5,5.5], [5.5,6], and [6,] in Å, respectively. Using these bins, the adjacency matrix has entries in an alphabet of size 14 (k=14), indicating bond type for covalently bonded atoms, and distance bin for all other atoms. We have trained the GG model on each target property individually. Further technical details are specified in the Supplementary materials sup .
III Results and discussion
III.1 Overview
We present an overview of the most relevant numerical results in Table 3. It contains the test errors for all combinations of regressors and representations and properties for models trained on 117 k molecules. The best models for the respective properties are (MG/GC), (MG/GG), (MG/GC), (MG/GC), (MG/GC), (BOB/KRR), (HDAD/KRR), (HDAD/KRR), (HDAD/KRR), and (BAML/RF). We do not show results for the other three energies, since identical observations as for can be made.
Overall, NN and KRR regressors perform well for most properties. The ML out-of-sample errors outperform DFT accuracy and reach chemical (target) accuracy, both defined alongside in table 3, for (GC), (BOB/KRR), (HDAD/KRR and MG/GG), (HDAD/KRR), and (BAML/KRR, MG/GC, HDAD/KRR, BOB/KRR, HD/KRR and MG/GG). For the remaining properties (, , , , and ) the best models come within a factor 2 of target accuracy, while all outperforming hybrid DFT accuracy.
In Fig. 2 out-of-sample errors as a function of training set size (learning curves) are shown for all properties and representations with the best corresponding regressor. It is important to note that all models on display systematically improve with training set size, exhibiting the typical linearly decaying behavior on a log-log plot Müller et al. (1996); Huang and von Lilienfeld (2016). Errors for most models shown decay with roughly the same slopes, indicating similar exponents in the power-law of error decay. Notable exceptions, i.e. property models with considerably steeper learning curves, are MG/GC for , MG/GG and HDAD/KRR for , CM/KRR and BoB/KRR for , HDAD/KRR and MG/GG for , and MG/GG for . These results suggest that the specified representations capture particularly well the effective dimensionality of the corresponding porperty in chemical space.
III.2 Regressors
Inspection of Table 3 indicates that the regressors can roughly be ordered by performance, independent of property and representation: GCGGKRRRFBREN. It is noteworthy how EN, BR, and RF regressors perform substantially worse than GC/GG/KRR. This can also be seen from the learning curves of all regressors presented in Figures S1 to S6 of the Supplementary Material. The performance of BR and EN improves only slightly with increased training set size and even gets worse for some property/representation combinations. These two regressors also exhibit very similar learning curves and BR preforms only slightly better than EN for most combinations. The only clear exception to this rule is for ZPVE and together with HDAD, where BR preforms significantly better than EN. Also, BR and EN errors rapidly converge to a constant w.r.t. training set size for all representations and properties, except for HDAD, which is the only representation which has a noteworthy improvement with increased training set size for some properties. The constant learning rates are not surprising as (a) the number of free regression parameters in BR and EN is relatively small and does not grow with training set size, and as (b) the underlying model is a linear combination with small flexibility. This behavior implies error convergence already for relatively small training sets.
RF performs poorly compared to GC, GG and KRR for all properties except for , the highest lying fundamental vibrational frequency in each molecule. For this property RF yields an astounding performance with out-of-sample errors as small as single digit cm*-1*. B3LYP achieves a mean absolute error of only 28 cm*-1* with respect to experiment Sinha et al. (2004). The distribution of , Figure 1 of reference Ramakrishnan and von Lilienfeld (2015), suggests a simple reason for this: There are three distinct peaks which correspond to typical C-H, N-H and O-H stretch vibrations in increasing order. Therefore the principal learning task in this property is to detect if there is an OH group, and if not if there is an NH group. If neither group is present, CH will yield the vibration with the highest frequency. As such, this is essentially about classifying which bonds are present in the molecule. RF works by fitting a decision tree to the target property. Each branch in the tree is based on an inequality of one entry in the representation. RF should therefore be able to identify which bonds are present in a molecule, simply by looking at the entries in the each element pair, and/or triplet bin of the representations. For RF, a fractional importance can be assigned to each input feature (the sum of all importances is 1.0). Analyzing the importance of the bins in HDAD of the RF model reveals that the three bins with highest importance are: O-H placed at 0.961 Å, N-H placed at at Å and C-C-H at radians with an importance of 0.587, 0.305 and 0.064 respectively. These three first bins constitute 96 of the prediction of and distances of the O-H and N-H bins are very similar to O-H and N-H bond lengths. C-C-H is placed on radians which means that it has to correspond to a linear C-C-H (alkyne) chain which implies that the two carbons must be bonded by a triple bond (typically the C-H with the lowest pKa and the highest C-H stretch vibration).
KRR performs remarkably well on average. For energetic properties, as well as , it yields the lowest overall errors in combination with HDAD and BoB, respectively. Its outstanding performance is not unsurprising in light of the multiple previous ML studies dealing with compositional as well as configurational spaces. The neural network flavors GC and GG, however, yield better performance on average, and the lowest errors for all electronic properties, i.e. dipolemoment, polarizability, HOMO/LUMO eigenvalues and gaps.
III.3 Representations
As one would expect, HDAD contains more relevant information and thus it always outperforms HD when using KRR. Tests also showed that an HDA representation systematically yields errors in between HDAD and HD, and similar observations hold for BR and EN regressor. In the case of RF, however, we observe little difference between HDAD and HD, and HD can even yield slightly lower errors than HDAD. In our opinion, this is due to the additional bins of angles and dihedrals rather adding noise than signal. By contrast, the separation of distances, angles and dihedral angles into different bins is not a problem for the KRR methods because the kernels used are purely distance based. This makes it possible for KRR to exploit the extra three- and four-body information in HDAD and to gain an advantage over HD. We note however that the remarkable performance of HDAD is possible despite its striking simplicity. As illustrated in Fig. 1 and discussed above, characteristic chemical behavior can be directly obtained by human inspection of HDAD. As such, HDAD corresponds to a representation very much ”Occam’s razor style”. Unfortunately, due to its discrete nature and its origin in sorting distances, HDAD will suffer from lack of differentiability, which might limit its applicability when modeling forces or other non-equilibrium properties.
MARAD, in principle containing the same information as HDA, performs similarly to BAML—yet, MARAD requires no prior knowledge about the physics encoded in the universal force-field. BoB and CM, while previously considered state of the art, result in relatively poor performance except for where BoB/KRR gives the lowest overall error. The simple bag of bonds, representing lists of all the pairwise Coulombic interactions appears to coincide ideally with the corresponding electronic spread in the molecule. Due to its inherent lack of uniqueness, however, we believe that also this property must converge to a finite error in the limit of infinite training set size. ECFP4 produces out-of-sample errors on par or slightly better than CM/KRR for intensive properties (, HOMO/LUMO eigenvalues and gap), however the model produces errors that are off-the-chart for all extensive properties (, , , and ).
IV Conclusions
We have benchmarked many combinations of regressors and representations on the same QM9 data set consisting of 131k organic molecules. For all properties, the best ML model prediction errors reach the accuracy of hybrid DFT with respect to experiment. For 8 out of 13 distinct properties (atomization energies, heat-capacity, , , ) out-of-sample errors reach levels on par with chemical accuracy, or better, using a training set size of 117k (90 of QM9 data set) molecules. For the remaining properties , , , , and , errors of the best models come within a factor 2 of chemical accuracy.
Regressors EN, BR, and RF lead to rather high out-of-sample errors, while KRR and graph based NN regressors compete for the lowest errors. We have found that GC, GG, and KRR have best performance across all properties, except for the highest vibrational frequency for which RF performs best. There is no single representation and regressor combination that works best for all properties (though forthcoming work with further improvements to the GG based models indicates best in class performance across all properties Gilmer et al. (2017)). For intensive electronic properties (, HOMO/LUMO eigenvalues and gap) we have found MG/NN to yield the highest predictive power, while HDAD/KRR corresponds to the most accurate model for extensive properties (, , , and ). The latter point is remarkable when considering the simplicity of KRR, being just a linear expansion of property in training set, and HDAD, being just histograms of distances, angels, dihedrals. Using BR and EN is not recommended if accuracy is of importance, both regressors perform worse across all properties. Apart from predicting highest fundamental vibrational frequency best, RF based models deliver rather unsatisfactory performance. The ECPF4 based models have shown poor general performance in comparison to all other representations studied; it is not recommended for investigations of molecular properties.
We should caution the reader that all our results refer to equilibrium structures of a set of only 131 k organic molecules. While 131k molecules might seem sufficiently large to be representative, this number is dwarfed in comparison to chemical space, i.e. the space populated by all theoretically stable molecules, estimated to exceed 1060 for medium sized organic molecules Kirkpatrick and Ellis (2004). Furthermore, ML models for predicting properties of molecules in non-equilibrium or strained configurations might require substantially more training data. This point is also of relevance because some of the highly accurate models described herewithin (MG based) currently use bond based graph connectivity in addition to distance, raising questions about the applicability to reactive processes.
In summary, for the organic molecules studied, we have collected numerical evidence which suggests that the out-of-sample error of ML is consistently better than estimated DFT accuracy While this does not prove that ML models would reach same error levels if more accurate, explicitly electron correlated or experimental reference data was used, there are studies which suggest that relevant property approximations do become smoother when using higher levels of theory Ghasemi et al. (2010); Ramakrishnan et al. (2015). As such, our study indicates that future reference data sets for training modern state-of-the-art machine learning models will require levels of theory which go beyond hybrid DFT. Finally, we note that future work could deal with improving representations and regressors with the goal of reaching similar predictive power using less data.
Acknowledgements.
The authors thank Dirk Bakowies for helpful comments, and Adrian Roitberg for pointing out an issue with the use of partial charges in the neural net models in an earlier version of this paper. O.A.v.L. acknowledges support from the Swiss National Science foundation (No. PP00P2_138932, 310030_160067), the research fund of the University of Basel, and from Google. This material is based upon work supported by the Air Force Office of Scientific Research, Air Force Material Command, USAF under Award No. FA9550-15-1-0026. This research was partly supported by the NCCR MARVEL, funded by the Swiss National Science Foundation. Some calculations were performed at sciCORE (http://scicore.unibas.ch/) scientific computing core facility at University of Basel.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136 , B 864 (1964).
- 2Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140 , A 1133 (1965).
- 3Burke (2012) K. Burke, J. Chem. Phys. 136 , 150901 (2012) . · doi ↗
- 4Koch and Holthausen (2002) W. Koch and M. C. Holthausen, A Chemist’s Guide to Density Functional Theory (Wiley-VCH, 2002).
- 5Cohen et al. (2012) A. J. Cohen, P. Mori-Sánchez, and W. Yang, Chem. Rev. 112 , 289 (2012) . · doi ↗
- 6Plata and Singleton (2015) R. E. Plata and D. A. Singleton, J. Am. Chem. Soc. 137 , 3811 (2015) . · doi ↗
- 7Medvedev et al. (2017) M. G. Medvedev, I. S. Bushmarinov, J. Sun, J. P. Perdew, and K. A. Lyssenko, Science 355 , 49 (2017) . · doi ↗
- 8Ramakrishnan et al. (2015) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, J. Chem. Theory Comput. 11 , 2087 (2015) . · doi ↗