Constant Size Molecular Descriptors For Use With Machine Learning

0.82 MB·

by  Christopher R. Collins

Constant Size Molecular Descriptors For Use
With Machine Learning
Christopher R. Collins,† Geoffrey J. Gordon,‡ O. Anatole von Lilienfeld,¶ and
David J. Yaron∗,†
7
†Department of Chemistry, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213
1
0 United States
2 ‡Machine Learning Department, Carnegie Mellon University, Pittsburgh, Pennsylvania

n 15213 United States
a
J ¶Department of Chemistry, Institute of Physical Chemistry and National Center for

3 Computational Design and Discovery of Novel Materials (MARVEL), University of Basel,
2
4056 Basel, Switzerland

]
h
p
-
m Abstract 1 Introduction
e
h
c A set of molecular descriptors whose length is Cheminformatics has a long tradition of using
.
s independent of molecular size is developed for the techniques of machine learning to predict
c
i machine learning models that target thermo- properties of molecules. A particular focus has
s
y dynamic and electronic properties of molecules. been properties that emerge from interactions
h These features are evaluated by monitoring per- of the molecule with complex external systems,
p
[ formance of kernel ridge regression models on such as arise in drug discovery.1,2 More recently,

well-studied data sets of small organic molecules. informatics has been successfully used to de-
1
v The features include connectivity counts, which velop models of the properties of the molecules
9
requireonlythebondingpatternofthemolecule, themselves. This includes thermodynamic prop-
4
6 and encoded distances, which summarize dis- erties,3–12 electronic13–26 properties and molecu-
6 tances between both bonded and non-bonded lar structure.27–30 Successful predictions of such
0
. atoms and so require the full molecular geome- molecular properties provides an alternative to
1
try. In addition to having constant size, these quantum chemical computations, leading to use-
0
7 features summarize information regarding the ful predictions at a small fraction of the cost of
1
local environment of atoms and bonds, such quantum computations.4,5 Because the target
:
v
that models can take advantage of similarities properties involve only the molecule itself, as
i
X
resulting from the presence of similar chemical opposed to interactions with complex external
r
fragments across molecules. Combining these systems, the descriptors that bring the most
a
two types of features leads to models whose per- explanatory power to the informatics tasks may
formance is comparable to or better than the differ from those developed for targets involving
current state of the art. The features introduced external systems. This paper explores a number
here have the advantage of leading to models of descriptors, or features, for predicting the
that may be trained on smaller molecules and properties of individual molecules.
then used successfully on larger molecules. A key dimension along which informatics mod-
els of molecular properties differ is the degree
to which results from chemical computations
are used as features for the prediction task. On
1
one end of this spectrum, the model corrects features, that count occurrences of patterns re-
predictions from standard computational chem- lated to three or more bonded atoms, are also
istry methods, so that they better match more explored.
expensive calculations17,18,31 or experiments.3,32 The second class of features are encoded dis-
At the other end of the spectrum, the only infor- tances. These summarize the pairwise distances
mation used is what one might get from a line between atoms in a molecule through a feature
drawing or SMILES33 string. At the latter level, vector whose length is independent of molecu-
there have been many examples of using simple lar size. The feature vector is a smoothed his-
counts of atoms and bonds to make predictions togram of distances between atoms, including
of molecular properties.5,34 The low computa- non-bonded atoms.
tional cost of such models makes them useful Both classes of features summarize informa-
for searching molecular databases.35 Their re- tion regarding the local environment of atoms
stricted applicability to potential energy minima and bonds in the molecule. The resulting fea-
or known bonding patterns; however, reduces turevectorsthereforecontaininformationonthe
their predictive power. fragments present in the molecule. In particular,
Between the above two extremes are models when two molecules possess the same molecular
that utilize the 3D molecular geometry to make fragment, the interactions amongst the atoms
predictions. Models that use features based on in that fragment make similar contributions to
the Coulomb matrix, such as the recent Bag of the feature vector. Use of these feature vectors
Bonds (BoB) feature vector, have been shown to predict similarity between two molecules, as
to be highly effective with this level of informa- in the kernel methods below, therefore incorpo-
tion.5,34 However, the length of these feature rates similarity resulting from the presence of
vectors scales formally as the square of the num- similar fragments. In addition, the feature vec-
berofatomsinthelargestmoleculeofinterest.36 tors include information regarding interactions
There have been attempts to rectify this issue, amongst the fragments. The features may there-
but they have been unable to achieve the same fore enable models to generalize from smaller to
level of performance.20 larger molecules.
A primary goal of the work presented here is
to develop feature vectors whose length is inde-
2 Related Work
pendent of molecular size and that subsequently
lead to models that generalize well from smaller
The features studied here are somewhat related
to larger systems. This should be possible if
tofeaturesdevelopedforneuralnetmodelsofpo-
we assume: a) an effective chemical periodicity
tentials.10,41–43 In those models, the total energy
where large systems consist of smaller fragments
is written as a sum over atomic energies, with
which resemble exemplary chemistries present in
the atomic energies being functions of the local
the data base used for training, and b) locality
environment.44 The environment is described us-
implying additivity. For certain systems, these
ing features based on two and three body corre-
assumptions can break down; in particular, the
lation functions around the atom. The encoded
latter one when strong electron delocalization
distance features explored here may be viewed
effects are present as in metals or in the π-bonds
as smoothed versions of two-body correlation
of conjugated polymers. To achieve our goal,
functions between each pair of elements. We
two classes of features are explored.
note that the latter lack uniqueness, a crucial
The first class of features are based on con-
feature property for guaranteeing error conver-
nectivity counts, which have been explored in
gencewithtrainingsetsize.20 Inparticular, they
a number of different ways.8,37–40 Here, these
will fail to distinguish homometric molecules.
features include counts of atom types and bond
Below, we use kernel-based learning methods
types, which we refer to as rank-1 and rank-2
to explore the utility of the features introduced
features since they summarize patterns related
here. Our studies therefore have some similari-
to one and two atoms, respectively. Higher-rank
2
ties to the Regularized Entropy Match Smooth atoms.46
Overlap of Atomic Positions (REMatch-SOAP) The QM7 data set corresponds to a subset
kernel, which provides an alternative means for of the GDB13 molecules, consisting of all 7,165
predicting molecular properties using only 3D molecules that contain 7 or fewer heavy atoms
molecular geometries.23 REMatch-SOAP does of elements C, N, O, or S for which optimized
not introduce a molecular feature vector and, structures and atomization energies were calcu-
instead, provides a means to directly compute lated with the Perdew-Burke-Ernzerhof hybrid
the similarity between two molecular structures, functional (PBE0).4 Models for the atomization
as needed for kernel methods. The approach be- energies in QM7 have been extensively stud-
gins by computing the similarity between each ied, primarily using features derived from the
of the atoms in the two molecules being com- Coulomb matrix.4–6,11,34 As such, this data set
pared, where the similarity measure includes the provides a useful benchmark for testing models
molecular environment. The molecular similar- and comparing to previous work.
ity is then obtained from a regularized entropy TheQM7bdatasetextendstheQM7dataset
match of the resulting atomic similarities. to molecules containing chlorine, and includes
In contrast to REMatch-SOAP, the features a broader range of properties such as HOMO
introduced here are not limited to kernel meth- and LUMO energies, excitation energies, and
ods and, when used in kernel methods, allow polarizabilities calculated at various levels of
us to separately consider the influences of three theory.16 QM7b allows performance of machine
factors: the feature vector, the distance metric learning models to be explored on a wider va-
used to connect features to the similarity be- riety of properties. Past works have also used
tween molecules, and the kernel function. Our QM7b to explore multi-target regression meth-
emphasis here is on the features themselves. We ods which, by simultaneously predicting various
use cross validation to select the kernel, includ- properties, can potentially benefit from relation-
ing the width parameter for an exponential ker- ships between the target properties.6,16,34
nel function and the choice of a Laplacian kernel, The QM9 data set47 contains data for a
which uses an L1 distance metric, versus a Gaus- subset of the GDB17 molecules, consisting of
sian kernel, which uses an L2 distance metric. 133,885 molecules that contain nine or fewer
The performance of the resulting models on well heavy atoms. (The name QM9 is used here to
studied data sets is examined, with particular beconsistentwiththatofQM7.) TheQM9data
attention paid to the ability of the method to set has been used in several studies19–21,34 and
scale between different size molecules. includes optimized structures and 18 different
properties calculated using B3LYP/6-31g(2df,p).
Below, QM9 is used to explore the degree to
3 Data Sets
which models training on smaller molecules can
transfer to larger molecules. The large amount
This work compares performance of various
of data in QM9 also allows more extensive stud-
molecular features on data sets that have been
ies on the degree to which inclusion of additional
the subject of past work, demonstrating the
data improves model performance.
predictive power of machine learning models
to predict molecular properties. These data
sets utilize molecules selected from the following 4 Machine learning
two data sets, each of which contains only the
bonding pattern as summarized in a SMILES A variety of machine learning methods have al-
string. The GDB13 data set contains 977 mil- ready been used to predict molecular properties,
lion molecules made of C, H, O, N, S, and Cl, including linear regression, linear ridge regres-
with up to 13 heavy atoms,45 and the GDB17 sion, support vector regression, kernel ridge re-
data set contains 166 billion molecules made gression, regression trees, k-nearest neighbors,
up of C, H, O, N, and F with up to 17 heavy and neural nets/deep learning.5,10,16 For this
3
study, we have elected to use two standard ma- The average of these MAEs is then taken as a
chine learning methods,48 linear ridge regres- measure of model performance.
sion (LRR) and kernel ridge regression (KRR) Here, model selection is based on a loop over
with the Gaussian and Laplacian kernels. These k=5 splits. Within this loop, one split is ig-
methods were chosen due to their relative sim- nored and the remaining 80% of data is used
plicity, compared for instance to deep learning, to evaluate each set of hyperparameters, via 5-
and because past work has shown them to be fold cross validation. This generates 5 figures
effective on the data sets studied here.5,34 of merit for each set of hyperparameters. The
The models involve very many parameters hyperparameters that have the best average per-
which, during “model training”, are adjusted formance are then selected as the final model.
to obtain predictions that agree with a given The MAE for this model, reported below and
set of examples. In addition, both LRR and in the Supporting Information, is the average
KRR involve hyperparameters that serve to de- from a final 5-fold cross validation performed on
fine further the model form. The wide range of the entire data set, using the selected hyperpa-
values for the hyperparameters leads to a wide rameters. (See Supporting Information for the
range of possible models. “Model selection” in- hyperparameters of each model.)
volves searching through these to find a model,
or equivalently a set of hyperparameters, that
5 Features
leads to good performance. Here, we do an ex-
haustive search over a grid of hyperparameter
The primary goal of this work is a compari-
values using cross validation as described below.
son of models developed using different types
InLRR,thereisonehyperparameter, α, which
of molecular features. Ideally, identical molecu-
is the weight on the ridge term. When α is zero,
lar geometries should lead to identical feature
the model produces the same result as linear
vectors, implying that the features are invari-
regression. Finite values for α penalize larger re-
ant with respect to molecular translations and
gressionparametersandsousingfinitevaluesfor
rotations, and to reordering of the atoms.
α favors simpler models (those with smaller re-
For models trained on smaller molecules to
gression parameters). For LRR, the following α
transfer to larger molecules, a number of other
values were scanned over: {10−5,10−3,...,101}.
properties of the feature vectors are desirable.
KRR retains the α parameter of LRR and adds
One such property is that the length of the fea-
two more. The choice of kernel for KRR can
turevectorbeindependentofmolecularsize. For
be considered as a discrete hyperparameter and
theCoulombmatrixandBoBfeatures,thisprop-
here the Gaussian (exp(−γ||x−y||2)) and Lapla-
2 erty does not hold, as the number of elements in
cian (exp(−γ||x−y|| )) kernels are included in
1
the feature vector scales quadratically with the
the search. For each of these kernels, there is an
number of atoms in the largest molecule in the
additional continuous hyperparameter, γ, which
data set. To support models of differently sized
sets the width of the kernel. The search over α
molecules, such feature vectors are expanded
and γ is exhaustive on a square grid with values
to a length that can accommodate the largest
{10−11,10−9,...,10−1}.
molecules in the data set, with zeros added to
Model selection and testing use k-fold cross
pad the extended regions for smaller molecules.
validation (CV). In k-fold CV, the examples
For the remaining features considered below,
in the data set are split into k equally sized
the feature vectors have a length that depends
sets. One of these folds is then held out, and
on the molecular diversity, e.g. atom and bond
the examples in the remaining folds are used to
types, but does not depend on molecular size.
train a model. The mean absolute error (MAE)
for predictions made on the examples in the
hold out set is then computed. This process
is repeated, with each one of the k folds being
the hold out set, leading to k different MAEs.
4
5.1 Null model catenated to make a single feature vector that
is invariant to reordering of the atoms. Because
A “null” model is used to provide a baseline
the BoB model is just a reordering of Coulomb
measure of the difficulty of the prediction task.
matrix elements, the length of the vector still
The null model always predicts the mean value
grows quadratically with the size of the largest
ofthetrainingdata. TheMAEofthenullmodel
molecule in the data set. BoB has also recently
therefore reflects the spread of the data.
been extended to include three-body and higher
terms,49 leading to feature vectors who length
5.2 Coulomb Matrix and Bag of growsasthethirdorhigherpowerofthenumber
Bonds of atoms in the largest molecule.
A class of features that has led to well-
5.3 Connectivity Counts
performing models of molecular properties
are those derived from the Coulomb ma-
This class of features summarizes the bonding
trix,4–6,16–19,21
pattern of a molecule through features that
(cid:40) count occurrences of some pattern regarding the
0.5Z2.4, if i = j
M = i (1) bonding between atoms. Each feature may be
ij ZiZj , if i (cid:54)= j assigned a rank, based on the number of distinct
|ri−rj|
atoms that are examined when testing for the
where Zi and ri are the atomic number and presence of the pattern. Rank-1 features count
Cartesian position of the ith atom. This feature atom types, such as element or element plus
is appealing in that it summarizes the molecular coordination number, leading to a vector whose
structure in terms of Coulomb interactions be- length is equal to the number of atom types.
tween atoms. Alternative forms for the interac- Rank-2 features count bond types (single, aro-
tion, other than Coulomb interactions between matic, double, triple), leading to a vector whose
bare nuclei, have been tried but were not found length is the number of bond types. Rank-3
to enhance model performance.34 features count triplets of bonded atoms, and so
While the Coulomb matrix is invariant to on.
molecular rotations and translations, it is not For rank 1, each element of the feature vector,
invariant to reordering of the atoms. A num- v , holds the number of atoms in the molecule
T
ber of means have been explored to address the that have type T,
dependence on atom order, including sorting
(cid:88)
the matrix based on the magnitude of the norm v = δ (2)
T A(i),T
of the rows, using randomly permuted sets of
i
matrices, and using the eigenvalues of the ma-
where the sum is over atoms, i, A(i) is a func-
trix.5,6,16 (The results generated here do not sort
tion that returns the atom type, and δ is the
the Coulomb matrix, and simply arrange the
Kronecker symbol. Two atom typing schemes
lower triangle of the symmetric matrix into a
are explored here, leading to two different rank-
vector.) A recent successful approach rearranges
1 feature vectors denoted as 1 and 1C. For
theCoulombmatrixintoa“BagofBonds”,BoB,
1, atom type is defined by element only. For
feature.34 In BoB, off-diagonal elements of the
1C, atom type is defined by a combination of
matrix are first divided into bags, based on the
element and coordination number. Coordina-
elements involved in the Coulomb interaction
tion number is taken as the number of bonds in
(CH, CC, CO, etc). The values in each bag are
which the atom participates.
then sorted from high to low values. The maxi-
For rank 2, each element of the feature vector,
mum size of each bag, across the molecules in
v , holds the number of bonds in the molecule
the data set, is then determined and zeros are T
added to each bag so that all molecules have
bags of the same length. The bags are then con-
5
that have type T, vector independent of the size of the molecules
included in the data set. We include only rank-2
(cid:88)
v = δ (3) encodings, which summarize information regard-
T B(i,j),T
i<j ing atom-atom distances. Extension to higher
ranks would lead to a substantial increase in the
where the sum is over atoms, and B(i,j) is a
length of the total feature vector and are not
function that returns the bond type between
explored here.
atoms i and j, or null if a bond is not present.
The encoded distance features may be viewed
Two bond typing schemes are explored here,
as generalizations of the above rank-2 connectiv-
leading to two different rank-2 feature vectors
ity features. Because the bond types of Eq. (3)
denotedas2and2B. For2,bondtypeisdefined
are based on atom-atom distances, the rank-2
bythetwoelementsparticipatinginthebond. A
connectivitycountsofSection5.3maybeviewed
bond is assumed to be present if the separation
as encoding atom-atom distances in a coarse
between atoms is less than the cutoffs listed in
grained manner. To introduce the notation we
the Supporting Information. For 2B, bond type
will use for encoding, we rewrite Eq. (3) as
is defined by the two elements participating in
the bond and the bond order, with bond orders (cid:88)
v = 1(r ∈ (rlow,rhigh]) (6)
assigned based on bond length as discussed in T jk T T
j<k
the Supporting Information.
For rank 3, each element of the feature vec- where rlow and rhigh are the lower and upper
T T
tor corresponds to a tuple of two bond types,
limitsofthedistancerangesthatdefinethebond
(T ,T ), with a value
1 2 type, T, for the corresponding pair of elements,
and r is the distance between atoms j and k.
(cid:88) jk
v = δ δ (4)
(T1,T2) B(i,j),T1 B(j,k),T2 Eq. (6) may be viewed as creating a histogram
i<j<k for each pair of elements. For feature vector 2,
the histogram has just one bin and reports the
where the sum is over atoms. The feature counts
number of bonds between those elements in the
only situations where atoms i and k are bonded
molecule. For feature vector 2B, the number of
to a common atom, j. The summation is thus
bins in the histogram is the number of allowed
over all triples of atoms to which a bond angle
bond orders between those elements.
would typically be ascribed. Higher-rank pat-
Asimpleapproachtoaddingmoreinformation
terns are generalizations of Eq. (4) to higher
to the feature vectors is to add more bins to the
order. For example, the rank 4 summation is
histogram. The information encoded can also
(cid:88) be expanded to include distances between non-
v = δ δ δ
(T1,T2,T3) B(i,j),T1 B(j,k),T2 B(k,l),T3 bonded atoms. To accomplish this, we create a
i<j<k<l
(5) uniformly spaced grid,
with the summation being over all sets of four
(rend −rstart)i
atoms to which a dihedral angle would typically d = rstart + (7)
i
N −1
be ascribed. grid
for i = 0,1,...,N −1 .
grid
5.4 Encoded distances
For each pair of elements, we then get a set of
The connectivity count features of Section 5.3 Ngrid features,
summarize the bonding pattern of a molecule
(cid:88)
but do not include information about the bond vi = 1(rjk ∈ (di,di+1]) . (8)
lengths, bond angles, or dihedral angles. The j<k
encoded features developed here allow informa-
The feature vector of Eq. (8) has two disadvan-
tion on the 3D structure of the molecule to be
tages. The first is that a small change in r
used while keeping the length of the feature jk
6
can lead to a discontinuous change in the fea- limitingthegeodesicdistance,G,betweenatoms
ture vector when r is near a grid point, d . included when computing the feature,
jk i
In addition, for values of N that are large
grid (cid:88)
compared to the number of pairs of elements v = f(r ,d ,β)1(G ≤ Gmax) (10)
i jk i jk
in the molecules of interest, the feature vector j<k
becomes sparse, making the learning prone to
where G is an integer specifying the number
overfitting. jk
of bonds along the shortest path between the
To overcome these disadvantages of the
atom j and atom k. The last term limits the
histogram-like features, we generalize Eq. (8)
summation to pairs of atoms that are separated
by replacing 1(r ∈ (d ,d ]) of Eq. (8) with
jk i i+1 by up to Gmax bonds. Other works have taken
an encoding function, f,
a similar approach by setting a cutoff based on
(cid:88) the Euclidean distance between atoms.10,11,41–43
v = f(r ,d ,β). (9)
i jk i
In the current studies, model performance was
j<k
found to be somewhat better when using a cut-
f returns a value in the range [0,1] and may off based on geodesic distance. Examples of
depend on a parameter, β, that is used to alter encoded features for benzene can be seen in
the smoothness of the function. We consider Figure 1.
two general classes of encoding functions: prob- The parameters specifying the encoded dis-
ability density functions (PDF) and cumulative tance features were tuned by evaluating perfor-
distribution functions (CDF). For both of these, mance on the QM7 data set. The Supporting
we consider three types of functions: the normal Information discusses the choice of N and β
grid
distribution, the logistic distribution, and the of Eq. (9). For the results reported below, β
spike distribution (Table 1). was 20 ˚A−1, and the grid consists of 100 points
between 0.2 ˚A and 6 ˚A. The choice of cutoff dis-
Table 1: Functions used for the encoded features tance, Gmax of Eq. (10), is based on the results
of Eq. (9). r is the distance between atoms, d
jk i in Figure 2, which show how model performance
are the grid points, β is a smoothing parameter, varies with Gmax for various encoded features.
x = β(r −d ), and (cid:15) is the distance between
jk i Results for the Coulomb matrix and BoB fea-
two grid points. tures are also shown, in which case Gmax is im-
plemented by setting to zero those values of the
Name Name2 f(r ,d ,β) in Eq. (9)
jk i Coulomb matrix that correspond to atoms sepa-
2NP Normal PDF exp(−x2)
rated by a geodesic distance greater than Gmax.
2NC Normal CDF 1+erf(x)
For nearly all features, the best performance is
2LP Logistic PDF exp(−x)/(1+exp(−x))2
obtained with Gmax of 2 or 3, corresponding to
2LC Logistic CDF 1/(1+exp(−x))
distances over which angles and dihedrals would
2SP Spike PDF 1(|r −d | < (cid:15)/2)
jk i be assigned to the molecular structure. This
2SC Spike CDF 1(r > d )
jk i is true of all features, when LRR is used, and
nearly all features, when KRR is used. The ex-
For the PDF class of functions, the feature
ceptions are KRR with 2SP (results not shown
vector may be viewed as adding noise to each
as the MAE is greater than 50 kcal/mol) and
of the atom-atom distances, r , with noise dis-
jk KRR with the Coulomb matrix. Note, however,
tributed according to the distribution of the
that when the Coulomb matrix elements are
encoding function. The CDF class of functions
reordered to form BoB features vectors, the op-
integrate the corresponding PDFs, leading to
timal Gmax is in the 2-3 range. The observation
a feature vector that increases monotonically
that Gmax is in this range for a broad range of
with distance and has a slope that indicates the
feature types and for both KRR and LRR mod-
presence of atom-atom separations near d .
i els suggests that the atomization energy has a
Empirical results, discussed below, suggest
local character, whereby the energy of a given
that model performance may be improved by
7
Figure 1: Encoded bond feature vectors for benzene. Left panels are 2NP and right panels are
2NC. Lower panels limit the geodesic distance between atoms, Gmax of Eq. (10), to 3. Upper panels
include all distances.
8
bond is primarily a function of its local envi- 4, with addition of rank 5 leading to a small
ronment. For the remainder of this work, Gmax degradation. Figure 3 also examines whether
is set to 3 for the encoded features. To allow inclusion of bond order improves performance,
comparison with past work, the BoB features e.g. whether 3B (red line) or 3 (black line) leads
do not limit Gmax. to a lower MAE. For feature vectors including
Because the KRR models outperform LRR by onlyconnectivitycounts, inclusionofbondorder
more than 0.5 kcal/mol, we also discuss only enhances performance at rank 3, but degrades
results from KRR. (Full results from both LRR performance at higher ranks. The lowest MAE
and KRR are tabulated in the Supporting Infor- for concatenated connectivity counts is 12B3B4
mation.) (3.40 kcal/mol).
The lower branch (dashed lines) of Figure 3
shows results from feature vectors that begin
6 Results
with 12LC, and so include encoded distances
from 3D molecular structures. For feature vec-
Throughout this section, we summarize results
tors that include encoded distances, inclusion
by listing the features, with notation as in Sec-
of rank 3 leads to substantial improvements in
tion 5, followed in parentheses by the MAE of
model performance and inclusion of rank 5 de-
a KRR model developed as in Section 4. For
grades performance. Inclusion of rank 4 either
example, 2B (7.69 kcal/mol) indicates that a
slightly improves or slightly degrades perfor-
KRR model using a feature vector 2B leads to
mance, depending on both the type of rank-2
anMAEof7.69kcal/mol. TheSupportingInfor-
encoding employed and the property being pre-
mation contains an extensive list of results from
dicted. In contrast to the connectivity only
both LRR and KRR models based on different
features, inclusion of bond order leads to better
feature sets. Selected results are discussed here.
performance at both rank 3 and rank 4. The
We first consider models that use a single
best performance on QM7 atomization energies
type of feature from Section 5 to predict the
is 12NP3B (1.19 kcal/mol). Inclusion of rank
atomization energies of the QM7 data set. For
4, as 12NP3B4B (1.25 kcal/mol), slightly de-
connectivity counts, the best performing single
grades performance. The performance is signifi-
features are 2B (7.69 kcal/mol) and 3B (6.36
cantly better than the 2.40 kcal/mol from BoB,
kcal/mol). The use of the 3D molecular struc-
which was the previous state-of-the-art for the
ture via encoded distance features lowers the
QM7 data set.34
error substantially, with the best feature being
We next use the QM9 data set to explore
2NP (2.45 kcal/mol). The other encoding func-
the extent to which models trained on smaller
tionsperformnearlyaswell, exceptfor2SP (146
moleculesmaybetransferredtolargermolecules.
kcal/mol)whichisessentiallyahistogramandso
KRR models, using the features of Table 2, were
gives a discontinuous and sparse feature vector.
trained to the 3,993 molecules of the QM9 data
This suggests that smooth encoding functions
set that have seven or fewer heavy atoms. Two-
lead to enhanced performance.
fold cross validation was used to select hyperpa-
Model performance can be enhanced by con-
rameters and the average MAE of the two folds,
catenating features of different ranks. The ef-
computed using the selected hyperparameters,
fects of adding features of increasing rank are
is listed in the column labeled CV (for cross val-
shown in Figure 3, with selected results shown
idation) in Table 2. The remaining columns list
in Table 2. The upper branch in Figure 3
the MAE as molecules from the QM9 data set
(solidlines)showsresultswhenonlyconnectivity
are added to the small molecules used to train
counts are included and all feature vectors begin
the model. Because the molecules are sorted
with 12B. These have the advantage of requir-
by size, the average molecular size increases as
ing only a line drawing or SMILES string of the
molecules are added. A useful quantitative mea-
molecule. Addition of higher-rank connectivity
sure of this increase is provided by the MAE of
counts improves performance up though rank
the null model, which reports on the spread of
9
Table 2: MAE of atomization energies for KRR models that use the listed feature vectors, for both
the QM7 and QM9 data sets. The QM9 model is trained to molecules with 7 or fewer heavy atoms
(3,993 molecules), yielding the cross validation error listed under CV. The remaining columns apply
the model to subsets of QM9 with increasing molecular size. The lowest MAE for each column is
shown in bold.
First X thousand
molecules of QM9
feature QM7 CV 20 50 133
Null 179.01 189.45 290.57 322.72 425.12
Coulomb matrix 3.37 3.83 43.78 66.56 106.09
BoB 2.40 2.43 9.68 22.21 29.99
1 14.58 15.46 18.20 21.65 20.99
2B 7.69 7.52 10.42 16.67 15.73
3B 6.36 6.84 8.26 15.39 15.73
12B 6.88 6.37 8.75 11.11 11.17
12B3 5.28 4.19 6.53 10.03 9.99
12B3B 4.51 4.25 6.34 8.36 8.59
12B34 3.64 3.68 5.70 9.18 9.89
12B34B 4.13 3.81 6.10 9.90 10.51
12B3B4 3.40 3.56 5.43 8.64 9.42
12B3B4B 3.87 3.69 5.92 9.53 10.25
2LC 2.74 2.10 3.26 5.83 5.76
2NP 2.45 1.65 2.65 5.18 4.71
12LC 1.68 1.57 2.42 4.20 4.09
12LC3 1.61 1.44 2.32 4.13 3.95
12LC3B 1.43 1.28 2.03 3.72 3.60
12LC34 1.73 1.38 2.21 4.36 4.21
12LC34B 1.46 1.34 1.89 3.77 3.52
12LC3B4 1.56 1.24 1.98 3.93 3.77
12LC3B4B 1.45 1.31 1.83 3.66 3.41
12NP 1.45 1.14 2.03 3.84 3.51
12NP3B 1.19 1.05 1.63 3.37 3.05
12NP3B4B 1.25 1.16 1.62 3.44 3.30
10