Chemical diversity in molecular orbital energy predictions with kernel ridge regression

TL;DR

This study evaluates kernel ridge regression for predicting molecular orbital energies across diverse datasets, highlighting the importance of molecular representation and demonstrating rapid predictions for new molecules.

Contribution

It compares two molecular representations in KRR and demonstrates the method's effectiveness across large, diverse datasets for predicting HOMO energies.

Findings

MBTR outperforms Coulomb matrix in accuracy

KRR achieves mean absolute error as low as 0.09 eV

Model successfully predicts energies for 10,000 unseen molecules

Abstract

Instant machine learning predictions of molecular properties are desirable for materials design, but the predictive power of the methodology is mainly tested on well-known benchmark datasets. Here, we investigate the performance of machine learning with kernel ridge regression (KRR) for the prediction of molecular orbital energies on three large datasets: the standard QM9 small organic molecules set, amino acid and dipeptide conformers, and organic crystal-forming molecules extracted from the Cambridge Structural Database. We focus on prediction of highest occupied molecular orbital (HOMO) energies, computed at density-functional level of theory. Two different representations that encode molecular structure are compared: the Coulomb matrix (CM) and the many-body tensor representation (MBTR). We find that KRR performance depends significantly on the chemistry of the underlying dataset and that the MBTR is superior to the CM, predicting HOMO energies with a mean absolute error as low as 0.09 eV. To demonstrate the power of our machine learning method, we apply our model to structures of 10k previously unseen molecules. We gain instant energy predictions that allow us to identify interesting molecules for future applications.

Figures

Tables

Table 1. Table 1 : Mean value, standard deviation and ranges of DFT-computed HOMO energies for all three datasets as depicted in Fig. 3 . QM9 and OE have similar mean values, ranges and standard deviations of HOMO energies, while the HOMO energies of the AA dataset spread out over a wide range of values due to metal cations with large atomic number that are attached to the amino acids and dipeptides.

Dataset	Mean [eV]	Std. dev. [eV]	Min [eV]	Max [eV]
QM9	-5.77	0.52	-10.45	-2.52
AA	-10.85	3.97	-19.63	-0.06
OE	-5.49	0.57	-12.70	-2.73

Table 2. Table 2 : Mean absolute errors (MAEs), root mean square errors (RMSEs) and R 2 coefficients for KRR predictions of molecules from the QM9, AA and OE datasets. Results are shown for the MBTR and CM representations. Errors and R 2 coefficients are measured on test sets of 10 k randomly selected out-of-sample molecules of each dataset, while 32 k molecules were used for training. We report errors and R 2 superscript 𝑅 2 R^{2} coefficients as the average over 10 repetitions, accompanied by the standard deviation of the mean. The corresponding reference-versus-predicted scatter plots are shown in Fig.  10 for the first run out of 10 repetitions.

	MAE [eV]		RMSE [eV]		R2
Dataset	MBTR	CM	MBTR	CM	MBTR	CM
QM9	0.086 $\pm 0.001$	0.151 $\pm 0.001$	0.118 $\pm 0.002$	0.207 $\pm 0.003$	0.950 $\pm 0.002$	0.845 $\pm 0.004$
AA	0.100 $\pm 0.001$	0.094 $\pm 0.002$	0.201 $\pm 0.006$	0.194 $\pm 0.007$	0.997 $\pm 0.001$	0.998 $\pm 0.001$
OE	0.173 $\pm 0.002$	0.336 $\pm 0.003$	0.239 $\pm 0.006$	0.435 $\pm 0.005$	0.821 $\pm 0.009$	0.413 $\pm 0.005$

Equations

$$C_{ij}=\{\begin{array}{ll}0.5Z_i^{2.4}& \text{if i=j}\\ \frac{Z_i{Z}_j}{\Vert {\mathbf{R}}_{\mathbf{i}}-{\mathbf{R}}_{\mathbf{j}}\Vert }& \text{if i=j}\end{array}.$$

$$k_{\text{Gaussian}}(\mathbf{x},{\mathbf{x}}^{\prime })=e^{-\frac{||\mathbf{x}-{\mathbf{x}}^{\prime }|{|}_2^2}{2{\sigma }^2}},$$

$$||\mathbf{x}-{\mathbf{x}}^{\prime }|{|}_2=\sqrt{\sum _i|{\mathbf{x}}_{\mathbf{i}}-{{\mathbf{x}}_{\mathbf{i}}}^{\prime }{|}^2}$$

$$k_{\text{Laplacian}}(\mathbf{x},{\mathbf{x}}^{\prime })=e^{-\frac{||\mathbf{x}-{\mathbf{x}}^{\prime }|{|}_1}{\sigma }}.$$

$$||\mathbf{x}-{\mathbf{x}}^{\prime }|{|}_1=\sum _i|{\mathbf{x}}_{\mathbf{i}}-{{\mathbf{x}}_{\mathbf{i}}}^{\prime }|.$$

$$\underset{\alpha }{\min }\sum _{i=1}^N(E^{\text{pred}}({\mathbf{x}}_i)-E_i^{\text{ref}}{)}^2+\lambda {\mathbf{\alpha }}^T\mathbf{K}\mathbf{\alpha },$$

$$\mathbf{\alpha }=(\mathbf{K}+\lambda \mathbf{I}{)}^{-1}{\mathbf{E}}^{\text{ref}}.$$

$$E^{\text{pred}}(\mathbf{x})=\sum _{i=1}^N{\alpha }_{i}k(\mathbf{x},{\mathbf{x}}_i).$$

Full text

I Introduction

Machine learning (ML) in molecular and materials science has gained increased attention in the last decade and its application domain is widening continuously.Rupp, von Lilienfeld, and Burke (2018) Applications include the search for improved and novel materials, Müller, Kusne, and Ramprasad (2016); Zunger (2018) computational drug design, (Ma et al., 2015) battery development, (Sendek et al., 2018; Shandiz and Gauvin, 2016) identification of new molecules for organic light-emitting diodes (Gómez-Bombarelli et al., 2016) or catalyst advancements for greenhouse gas conversion. (Goldsmith et al., 2018; Meyer et al., 2018) In the context of ab initio molecular science, ML has been applied to learn a variety of molecular properties such as atomization energies, (Hansen et al., 2013; Rupp et al., 2012; Huang and von Lilienfeld, 2016; Faber et al., 2017, 2018; Bartók et al., 2017; Collins et al., 2018; De et al., 2016; Ramakrishnan et al., 2015a; Montavon et al., 2013) polarizabilities, (Faber et al., 2018; Schütt et al., 2018; Huang and von Lilienfeld, 2016; Collins et al., 2018; Montavon et al., 2013) electron ionization energies and affinities, (Huang and von Lilienfeld, 2016; Collins et al., 2018; Montavon et al., 2013; Ghosh et al., ) dipole moments, (Pereira and de Sousa, 2018; Faber et al., 2017, 2018; Schütt et al., 2018; Huang and von Lilienfeld, 2016; Bereau, Andrienko, and von Lilienfeld, 2015; Bereau et al., 2018) enthalpies, (Faber et al., 2017; Huang and von Lilienfeld, 2016; Ramakrishnan et al., 2015a) band gaps, (Faber et al., 2017, 2018; Schütt et al., 2018; Huang and von Lilienfeld, 2016), binding energies on surfaces Todorović et al. (2019) as well as heat capacities. (Faber et al., 2017, 2018; Huang and von Lilienfeld, 2016) A few studies addressed the prediction of spectroscopically relevant observables, such as electronic excitations, (Huang and von Lilienfeld, 2016; Collins et al., 2018; Montavon et al., 2013; Pronobis et al., 2018) ionization potentials, (Huang and von Lilienfeld, 2016; Collins et al., 2018; Montavon et al., 2013; Ramakrishnan et al., 2015b) nuclear chemical shifts, (Rupp, Ramakrishnan, and von Lilienfeld, 2015) atomic core level excitations (Rupp, Ramakrishnan, and von Lilienfeld, 2015) or forces on atoms. (Rupp, Ramakrishnan, and von Lilienfeld, 2015) In this paper, we employ kernel ridge regression (KRR) to predict the energy of the highest occupied molecular orbital (HOMO), which is of particular current interest for the development of new substances and materials. Electronic devices based on organic compounds are widely used in the technological industry and frontier orbital energies give important information about the optoelectronic properties of possible candidate materials.

While the majority of ML studies focuses on ground state properties, a handful of studies have addressed the pre-

diction frontier orbital energies, for example using neural networks, (Montavon et al., 2013; Pyzer-Knapp, Li, and Aspuru-Guzik, ; Schütt et al., 2018; Faber et al., 2017; Ghosh et al., ) random forest models, (Pereira et al., 2017) or kernel ridge regression. (Huang and von Lilienfeld, 2016; Collins et al., 2018; Faber et al., 2018; De et al., 2016) To our knowledge, the best prediction accuracy reported for HOMO energy predictions on the well-known QM9 dataset (Ramakrishnan et al., 2014) of small organic molecules was achieved with a deep neural network and featured a mean absolute error (MAE) of 0.041 eV. (Schütt et al., 2018) QM9 was also employed by many other studies to explore the effect of molecular descriptors on prediction accuracy. (Rupp et al., 2012; Huang and von Lilienfeld, 2016; Faber et al., 2017; Schütt et al., 2018)

While these results are quantitative and valuable, it is not clear to what extent the predictive power of the employed methodology is transferable to other molecular datasets. Motivated by optoelectronic applications, we are interested in HOMO predictions for large optically-active molecules with complex aromatic backbones and diverse functional groups which differ notably from the QM9 dataset.

We here employ KRR on three different datasets, two of which have not yet been used for the prediction of molecular orbital energies with ML. These two less-known datasets consist of 44 k conformers of proteinogenic amino acids from a public database of oligo-peptide structures (Ropo et al., 2016) and 64 k organic molecules extracted from the organic crystals of the Cambridge Structural Database.(Schober, Reuter, and Oberhofer, 2016; Groom et al., 2016) In addition, we also use the well-known QM9 benchmark database of 134 k small organic molecules (Ruddigkeit et al., 2012; Ramakrishnan et al., 2014) as a third dataset to compare results with previous studies in this field on a common basis. For all three datasets, we have calculated reference HOMO energies with density functional theory (DFT).

Moreover, we compare the performance of two different molecular representations: the well-studied Coulomb matrix (CM), which is simple, easy to compute and yields fast and inexpensive ML predictions. The second is a constant-size representation recently introduced by one of us (Huo and Rupp, 2017) that relies on interatomic many-body functions including bonding and angular terms, called the many-body tensor representation (MBTR). While previous studies already demonstrated that the CM can easily be outperformed by more sophisticated molecular descriptors Huang and von Lilienfeld (2016); Collins et al. (2018); Faber et al. (2017), we aim to analyze the degree of accuracy that can be achieved with the simple and cheap CM in comparison to the costlier MBTR.

The primary goal of our study is the comparison of KRR performance across three datasets with different chemical diversity. We show that the accuracy of HOMO energy predictions with KRR depends – besides the choice of molecular representation – crucially on the chemistry of the underlying dataset. Our measurable acceptable accuracy for HOMO energy predictions is 0.1 eV. Experiments for HOMO energy determination typically have a resolution of several tenth of eV, and prediction errors of state-of-the-art theoretical spectroscopy methods commonly range between 0.1 and 0.3 eV. We demonstrate how differences in model performance across chemically diverse settings can be related to certain dataset properties. Moreover, we quantify molecular orbital energy predictions that are presently available for realistic datasets of technological relevance.

Once trained, our KRR model can make instant HOMO energy predictions for numerous unknown molecules at no further cost. We demonstrate this by producing a spread of HOMO energy predictions for a new dataset of 10k organic molecules Ramakrishnan et al. (2015a), whose original HOMO energies are unknown. With instant energy predictions for all 10k molecules we gain a rough estimate of the HOMO energy distribution for this dataset. We can further identify interesting molecules within a certain energy range for additional analysis. Hence, large numbers of new molecules, whose orbital energies have not yet been measured or computed, can be quickly screened for their usability in future applications. In this way, KRR can complement conventional theoretical and experimental methods to greatly accelerate the analysis of materials.

The manuscript is organized as follows: In Section II we introduce the three datasets used in this work. In Section III we briefly review the two descriptors that we employ in our ML framework, which is described in Section IV. In Section V we present our results and then discuss our findings in Section VI.

II Datasets

We train and evaluate our model on three different datasets, for which example molecules are depicted in Fig. 1. For all three datasets, we performed DFT calculations with the FHI-aims code. (Blum et al., 2009; Havu et al., 2009; Levchenko et al., 2015; Ren et al., 2012) We optimized the atomic structure of all molecules using the Perdew-Burke-Ernzerhof (PBE) functional (Perdew, Burke, and Ernzerhof, 1996) including Tkatchenko-Scheffler van der Waals corrections (PBE+vdW), (Tkatchenko and Scheffler, 2009) tight computational settings and the tier 2 basis sets of FHI-aims. For reasons of computational tractability, we also calculate HOMO energies with DFT, by taking the eigenvalue of the highest occupied molecular state from the PBE+vdW calculation.

Although DFT Kohn-Sham energies have limited accuracy, they provided us with a large, convenient and consistent dataset for developing our methodology. In the future, we will extend our study to HOMO energies computed with the more appropriate $GW$ method.Hedin (1965); Rinke et al. (2005) However, at present it is not possible to calculate hundreds of thousands of molecules with the $GW$ method with reasonable computational resources. In the following, we describe the datasets in more detail.

II.1 QM9: 134 k small organic compounds

This dataset is extracted from the QM9 database and consists of 133,814 small organic molecules with up to 9 heavy atoms made up of C, N, O and F atoms (Ramakrishnan et al., 2014). It contains small amino acids and nucleobases, pharmaceutically relevant organic building blocks, for a total of 621 stoichiometries. The QM9 database has been used in a variety of ML studies (Collins et al., 2018; Huang and von Lilienfeld, 2016; Bartók et al., 2017; Ramakrishnan et al., 2015a; Ramakrishnan and Lilienfeld, 2017; Schütt et al., 2018; Faber et al., 2017; Pronobis et al., 2018; Schütt et al., 2017; Faber et al., 2018; Lubbers, Smith, and Barros, 2018; Unke and Meuwly, 2018) and has become the drosophila of ML in chemistry.

II.2 AA: 44 k amino acids and dipeptides

This dataset, denoted AA, contains 44,004 isolated and cation-coordinated conformers of 20 proteinogenic amino acids and their amino-methylated and acetylated (capped) dipeptides. (Ropo et al., 2016) The molecular structures are made of up to 39 atoms including H, C, N, O, S, Ca, Sr, Cd, Ba, Hg and Pb. The amino acid conformers reveal different protonation states of the backbone and the sidechains. Furthermore, amino acids and dipeptides with divalent cations (Ca2+, Ba2+, Sr2+, Cd2+, Pb2+, and Hg2+) are included. Since all amino acids share a common backbone the complexity of this dataset lies in differing sidechains and differing dihedral angles. AA has been used to benchmark several ML models (Bartók et al., 2017; De et al., 2016; Artrith, Urban, and Ceder, 2017) and clustering techniques. (De et al., 2017)

II.3 OE: 64 k opto-electronically active molecules

This dataset, referred to as OE, consists of 64,710 large organic molecules with up to 174 atoms extracted from organic crystals in the Cambridge Structural Database (CSD). (Groom et al., 2016) Schober et al. have screened the CSD for monomolecular organic crystals with the objective to identify organic semiconductors with high charge carrier mobility. Schober, Reuter, and Oberhofer (2016); Kunkel et al. (2019a, b) For this study, we extracted molecules from the crystals and relaxed them in vacuum with the aforementioned computational parameters. The OE dataset is not yet publicly available. OE offers the largest chemical diversity among the sets in this work both in terms of size as well as number of different elements (Fig. 2). It contains the 16 different element types H, Li, B, C, N, O, F, Si, P, S, Cl, As, Se, Br, Te and I. The electronic structures are more complex than in QM9 and AA, containing, e.g. large conjugated systems and unusual functional groups.

II.4 Comparison of datasets

Figure 2 illustrates the chemical diversity present in our datasets. The molecular size distribution of the three datasets in Fig. 2 shows that QM9 and AA both contain molecules of similar sizes. The QM9 distribution exhibits a distinct peak at around 18 atoms, whereas AA has a bimodal distribution centered around 27 and 35 atoms. Conversely, the size distribution of OE is much broader, extending to molecules with as many as 174 atoms. In terms of element diversity, all datasets overlap on the 4 elements of organic chemistry: H, C, O, and N, as illustrated in Fig. 2. QM9 contains only F in addition, whereas AA branches out into common cations. OE offers the largest element diversity, including common semiconductor elements.

To compare our target property across the different datasets, we show distributions of the DFT pre-computed HOMO energies for each dataset in Figure 3. HOMO energies in QM9 and OE are centered around -6 eV, with a standard deviation of 0.5 eV and 0.6 eV, respectively. For AA, only a fraction of the HOMO energies are centered around -6 eV, while most HOMO energies are distributed over a wider range between -19.6 eV and -8 eV. The HOMO energies around -6 eV correspond to amino acids and dipeptides in bare organic configurations (free of cations). The HOMO energies between -19 and -8 eV correspond to amino acids and dipeptides with one of the six cations Ca2+, Ba2+, Sr2+, Cd2+, Pb2+ or Hg2+. The metal ions in AA shift the HOMO energies of the amino acids towards lower values.

Judging by the distributions of molecular size and element types in Fig. 2 and by the distributions of HOMO energies in Fig. 3, we expect QM9 to be learned relatively easily. QM9 is clustered both in chemical space (in terms of molecular size and element types) and in target space (in terms of HOMO energies). We therefore expect KRR to benefit from mapping similar input structures to similar target values. For OE, the target energy distribution is similar to QM9, but the molecular structures are widely spread through chemical space. Similarity is therefore present only in one space, while the other is diverse. It might be challening for KRR to map from the diverse chemical space of OE to its confined target space and we expect learning to be slow. The AA set is spread out both in chemical space and target space. As it turns out, this will not be a problem for learning.

In Fig. 4 the three datasets are visualized with the t-distributed stochastic neighbor embedding (t-SNE) dimensionality reduction technique for both the their CM and MBTR molecular representations. (van der Maaten and Hinton, 2008) From each dataset, we randomly picked 3,000 molecules, which are mapped to a two dimensional space by t-SNE. Based on the pairwise similarities between descriptors, we can identify patterns within the datasets. In Fig. 4a), we can see that the AA and QM9 sets form distinct clusters and that the two clusters (orange and blue) have almost no overlap in the CM representation. The MBTR t-SNE analysis produces the same result, as indicated by Fig. 4b) (cluster orientation is arbitrary). This implies that the molecules in QM9 are very different from those in AA and that we increase the chemical diversity of our study by including the AA dataset. Fig. 4 further illustrates that the OE set is maximally diverse in itself. The corresponding red point cloud is evenly spread over the whole figure. In contrast to the QM9 and AA sets, OE offers a variety of rigid backbones (also denoted scaffolds) to which different functional groups are attached. Chemical diversity in part arises due to the rich combinatorial space emerging from scaffold-functional group pairings. For example, OE contains molecules with conjugated and aromatic backbones and electron withdrawing and donating functional groups that are of technological relevance and completely absent from QM9 and AA (Kunkel et al., 2019a, b).

To quantify the similarity between molecules in our three datasets, we computed Euclidean distances for molecular pairs – represented by the CM – for 3,000 randomly chosen molecules for each dataset, as shown in Fig. 5. Molecular distances within QM9 are small, indicating great similarity among QM9 molecules. The OE distances are distributed evenly over a larger and wider range of distances, indicating high dissimilarity among OE molecules. The AA distance distribution has several separated clusters up to very large distances. The distances are separated and ordered by the atomic number of the cations, where each cluster corresponds to amino acids with different metal ions attached. The first cluster includes distances from bare amino acids to amino acids with Ca2+, followed by clusters including distances to amino acids with Sr2+, Cd2+, Ba2+, Hg2+ and Pb2+. Structural dissimilarity in the AA dataset mainly arises due to amino acids with different metal ions, while amino acids within the same cluster are highly similar to each other.

III Molecular representation

For the ML model to make accurate predictions, it is important to represent the molecules for the machine in an appropriate way. (Rupp, 2015; Bartók, Kondor, and Csányi, 2013; von Lilienfeld et al., 2015; Huo and Rupp, 2017) Cartesian (x,y,z)-coordinates, which are, for example, used for DFT calculations, are not applicable, since they are not invariant to translations, rotations, and reordering of atoms. In this work, we compare the performance of two different molecular representations.

III.1 Coulomb matrix (CM)

In the CM formalism, (Rupp et al., 2012) each molecule is represented by a matrix C,

[TABLE]

The CM encodes nuclear charges $Z_{i}$ and corresponding Cartesian coordinates $R_{i}$, with off-diagonal elements representing Coulomb repulsion between atom pairs and diagonal elements encoding a polynomial fit of free-atom energies to $Z$. An example of the CM for a molecule of OE is shown in Fig. 6. To enforce permutational invariance, we simultaneously sort rows and columns of all CMs with respect to their $\ell^{2}$-norm.

III.2 Many-body tensor representation (MBTR)

The MBTR (Huo and Rupp, 2017) can be viewed as a many-body expansion of the Bag of Bonds (BoB) representation, (Hansen et al., 2015) which in turn is based on elements of the CM. One-body terms of the MBTR describe the atom types that are present in the molecule. Two-body terms encode inverse distances between all pairs of atoms (bonded and non-bonded), separately for each combination of atom types. The inverse distances are sorted by increasing order and broadened into a continuous Gaussian distribution of inverse distances as shown in Fig. 7a). Three-body terms encode angle distributions for any triplets of atoms present in the molecule, as shown in Fig. 7b). Each $N$-body term has a broadening parameter (in total $\rho_{1}$, $\rho_{2}$ and $\rho_{3}$) that controls the smearing of atom type distribution, inverse distance distribution and angle distribution, respectively, and need to be fine-tuned for optimal KRR performance. We use the DScribe package dsc to compute the MBTR and the qmmlpack package to refine the MBTR hyperparameters for small training set sizes of up to 4k molecules. Exponential weighting was employed for the computation of inverse distance and angle terms. Here, we apply only two-body and three-body terms in the MBTR, since we found that the inclusion of one-body terms does not improve the performance, but increases computational time (we refer to Appendix B for details).

IV Machine learning method

IV.1 Kernel ridge regression

We employ kernel ridge regression (Hastie, Tibshirani, and Friedman, 2009) (KRR) to model the relationship between molecular structures and HOMO energies. In KRR, training samples are mapped into a high-dimensional space using nonlinear mapping, and the structure-HOMO relationship is learned in the high-dimensional space. This learning procedure is conducted implicitly by defining a kernel function, which measures the similarity of training samples in the high-dimensional space by employing a kernel function. In this work, we use two different kernel functions. The first kernel is the Gaussian kernel

[TABLE]

where

[TABLE]

is the Euclidean distance and $\bm{x}$, $\bm{x}^{\prime}$ are two training molecules represented by either the CM or the MBTR. We note that the Euclidean distance distribution between molecular pairs shown in Fig. 5 gives us a direct insight into the learning process of a Gaussian kernel function. The second kernel is the Laplacian kernel

[TABLE]

which uses the 1-norm as similarity measure,

[TABLE]

In Eqs. (2) and (4), $\sigma$ is the kernel width.

In the KRR training phase with $N$ training molecules, the goal is to find a vector $\bm{\alpha}\in\mathbb{R}^{N}$ of regression weights $\alpha_{i}$ that solves the minimization problem

[TABLE]

where the analytic solution for $\bm{\alpha}$ is given by

[TABLE]

The matrix $\mathbf{K}$ is the kernel matrix, whose elements represent inner products between training samples in the high-dimensional space, calculated as kernel evaluations $K_{i,j}:=k(\bm{x_{i}},\bm{x_{j}}$). The scalar $\lambda$ is the regularization parameter, which penalizes complex models with large regression weights over simpler models with small regression weights. The expression $E^{\textnormal{ref}}$ denotes the reference HOMO energy computed by DFT and $E^{\textnormal{pred}}$ is the predicted HOMO energy, which is obtained as sum over weighted kernel functions

[TABLE]

The sum runs over all molecules $\bm{x}_{i}$ in the training set with their corresponding regression weights $\alpha_{i}$. After training, predictions are made for test molecules that were not used to train the model, employing eq. (8) to predict out-of-sample molecules and estimate performance of the model. No further scaling or normalization of the data was done to preserve the meaning of the HOMO energies.

IV.2 KRR training, cross-validation and error evaluation

For each dataset QM9, AA and OE, we randomly selected a subset of 32 k molecules for training and a further 10 k molecules for out-of-sample testing. In order to obtain statistically meaningful results for the training and testing performance of KRR, we repeated the random selection of training set and test set nine more times for each dataset after the data was reshuffled. As a result, we acquired 10 different training sets of 32k molecules and 10 different test sets of 10k molecules for each dataset.

From each training set of 32k, we randomly drew 6 different subsets of sizes 1k, 2k, 4k, 8k and 16k, where smaller sets were always subsets of larger ones. For each of these subsets, we trained and cross-validated a KRR model, as described in the following paragraph. We then evaluated the KRR model on the corresponding test set of 10k by predicting HOMO energies for 10k out-of-sample molecules and by computing the MAE $\frac{1}{n}\sum_{i=1}^{n}|E^{\textnormal{pred}}_{i}-E^{\textnormal{ref}}_{i}|$ between predictions and DFT reference energies. For each dataset and for each training set size, we computed the average MAE value and its standard deviation across the 10 randomly drawn training and test sets. We then plotted all of these average MAEs as a function of training set size and as a result, attain one learning curve for each dataset, as shown in Fig. 8 and Fig. 11. We note that for the two molecular representations, CM and MBTR, identical training and test sets were used.

In the scope of this study, there are 3 types of hyperparameters:

1. (i)

the kernel function (Gaussian type or Laplacian type) 2. (ii)

MBTR hyperparameters (broadening values $\rho_{1}$, $\rho_{2}$ and $\rho_{3}$) 3. (iii)

KRR hyperparameters (kernel width $\sigma$ and regularization parameter $\lambda$)

The large number of possible values for the hyperparameters leads to a wide range of possible KRR models. The choice of kernel function depends on the molecular descriptor. The Gaussian kernel performs best on the MBTR, whose $N$-body terms themselves consist of Gaussian distributions. The Laplacian kernel, on the other hand, can better model piecewise smooth functions, such as discontinuities of the sorted CM. Therefore, we chose the Laplacian kernel for the CM and the Gaussian kernel for the MBTR. We provide learning curves in Appendix B to prove that we picked the optimal kernel for each descriptor.

The MBTR and KRR hyperparameters were simultaneously optimized in a cross-validated grid search for each of the three smaller training sets of 1k, 2k and 4k. In particular, MBTR hyperparameters $\rho_{1}$, $\rho_{2}$ and $\rho_{3}$ were varied on a grid of 10 points between $10^{-4}$ and $10^{-1}$ and KRR hyperparameters $\sigma$ and $\lambda$ were varied on a grid of 12 points between $10^{-12}$ and $10^{-1}$. The aim was to find the combination of best KRR and MBTR parameters. For more details on the cross-validated grid search we refer to Appendix C. The optimized model was then evaluated on 10k out-of-sample molecules from the test set. Finally, the MAE on the test set was reported as a data point on the learning curve for a given training set size.

We found that the values of the optimized MBTR hyperparameters do not change throughout the three small training sets. We therefore fixed the optimized MBTR hyperparameters for the larger larger training sets of 8k, 16k and 32k and optimized only the KRR hyperparameters (on a grid from $10^{-12}$ to $10^{-1}$).

Another measurent for the performance of a ML model is the $R^{2}$ coefficient, which describes the proportion of variability in a dataset that can be explained by the model. Although we do not use the R2 coefficient to optimize KRR hyperparameters, we employ this metric to interpret how well our models fit the data.

V Results

Upper panels in Fig. 8 show out-of-sample MAEs as a function of training set size ("learning curves”) for the different datasets and for CM and MBTR as descriptors. As expected, the MAE decreases for all datasets with increasing training set size (see also Fig. 11). The best MAEs for a training set size of 32k molecules are presented in Table 2. The lowest MAE is achieved for QM9, closely followed by AA. In constrast, the MAEs for OE are approximately twice as high. The learning rate (slope of the MAE curves) is highest for AA and lowest for OE, and is independent of the descriptor. The MBTR performs significantly better than the CM for QM9 and OE, while for AA, the learning curves of the two descriptors are the same within statistical errors.

The lower panels of Fig. 8 show the squared correlation coefficient $R^{2}$ as a function of training set size. $R^{2}$ increases systematically with training set size, but its rate varies across the three datasets. For AA it is close to 1 already for small training set sizes, whereas for QM9, $R^{2}$ starts off low and then approaches 1 for 32 k. For the OE dataset, $R^{2}$ steadily increases and reaches a value of 0.81 at 32 k.

Correlating $R^{2}$ and MAE, we observe that model fitting for OE is consistently poor (low $R^{2}$) and the prediction errors are consistently high (high MAEs). The CM appears to be less suitable for this prediction task.

Fig. 10 presents scatter plots for a fixed training set size of 32 k, which describes how well the KRR predictions correlate with the reference values in the test set. We observe the best correlation (i.e. predictive power) for AA and the worst for OE. The MBTR appears to have higher predictive power than CM for QM9 and OE, whereas MBTR and CM perform similarly for AA.

Fig. 11 summarizes the learning curves for the three datasets. The top panel shows the MBTR and the bottom panel the CM results. The MBTR-based KRR models generally produce faster learning rates and lower MAEs. The AA learning rate is particularly fast, and eventually leads to MAEs of 0.10 eV for MBTR and 0.09 eV for CM, which are comparable to those of QM9. The prediction quality for the OE dataset is notably worse in relation to the other datasets.

VI Discussion

VI.1 Dependence of KRR performance on dataset diversity

The learning success of KRR depends on the structural complexity of individual molecules (e.g. number of atoms, atom types, backbone types etc.) as well as on the diversity and redundancy within a dataset. Redundancy usually means that certain structural features occur frequently in a dataset, i.e. many data points are similar to each other. Diversity implies the opposite: Few instances in a dataset are similar to the rest. Redundant datasets are learned well with ML, even when trained on small portions. Diverse datasets can pose a problem for ML, even when applied to large data.

The differences in the learning curves we observe in Figs. 8 and 11 reflect the chemical differences in the three datasets employed in our study. QM9 is greatly redundant and includes molecules with simple structures. Therefore, it can be learned well even on small training set sizes, where redundancy is low.

The AA dataset has inbuilt redundancy, but also includes many different metal cations. For small training set sizes, where not enough similar structures per metal cation are present, the error is high. This situation then improves quickly with increasing traing set size. As a result, the learning rate of AA is faster than of QM9.

OE is highly diverse and includes molecules with complex structures. It has similar chemical element range as AA with 16 different element types, but features larger molecules and more structural diversity, as discussed in Section II.4. The diversity explains the high errors throughout all training set sizes and the slower learning rate. The t-SNE analysis in Fig. 4 confirms that OE has little overlap with QM9, regardless of the molecular representation. OE is the most diverse of the three datasets, since the corresponding point cloud does not cluster in a particular region and instead fills the whole space of Fig. 4.

The KRR learning rate is slowest for OE, due to the aforementioned chemical and structural complexity in this dataset. At a training set size of 32 k, the MAE is still twice as high as for QM9 and AA. With 0.173 eV, the MAE for OE is still too high for the spectroscopic applications we intend, which typically require errors below 0.1 eV. Advanced machine learning methodology and more sophisticated materials descriptors both help to reduce overall prediction errors (Montavon et al., 2012; Faber et al., 2017, 2018; Schütt et al., 2018; Pyzer-Knapp, Li, and Aspuru-Guzik, ) and we may test alternative approaches in further work.

Lastly, our study of three different datasets illustrates that MAEs are not transferable across datasets, even if they are evaluated for the same machine learning method and the same descriptor. If we had based our predictive power expectations for the OE set on the KRR QM9 performance, we would have been disappointed to find much larger errors in reality. It is therefore paramount to further investigate the performance of machine learning methods across chemical space.

VI.2 Dependence of KRR performance on molecular descriptor

Next, we discuss the relative performance for the CM and MBTR molecular descriptors. Overall, MBTR outperforms the CM across the datasets, which is in line with previous findings. (Huang and von Lilienfeld, 2016; Faber et al., 2017; Pereira et al., 2017; Collins et al., 2018) This is reasonable in light of the higher information content about atom types, their bond lenghts and angles encoded in the MBTR when compared to the CM matrix. The CM and the MBTR exhibit the same performance only for the AA dataset of molecular conformers, where complexity is dominated by the torsional angles, and bonding patterns are similar. This result is partly explained by the exclusion of torsional angle information into MBTR (four-body terms), while consistent chemical information in AA benefits the performance of the CM.

It is interesting that CM generally produces MAEs only twice as large as the MBTR with much smaller data structures and at a fraction of the computational cost. The CM representation is simple to compute, supplies benchmark results comparable to previous work and may prove a convenient tool for preliminary studies of large unknown datasets.

VI.3 Application of QM9 model

Our results for the QM9 dataset allow us to compare our findings with previous studies. Given a KRR training set of HOMO energies for 32 k molecules, we obtained MAE values of 0.086 eV and 0.151 eV with the MBTR and CM descriptors respectively. In 2017, Faber et al. (Faber et al., 2017) reported KRR results on 118 k QM9 molecules, using a molecular descriptor based on interatomic many-body expansions including bonding, angular and higher-order terms, (Huang and von Lilienfeld, 2016) which is comparable to the MBTR used in this work. The HOMO energy was predicted with an out-of-sample MAE of 0.095 eV, and the CM representation achieved an MAE of 0.133 eV. Our QM9 results are in very good agreement with this study, even if our KRR training set is much smaller. These errors are relatively small and comparable to experimental and computational errors in HOMO determinaton, which typically range in between several tenth of eV. Errors in HOMO energy predictions with machine learning may be further reduced by developing customized deep learning neural network architectures, which have been reported to produce an MAE of 0.041 eV after training on 110 k QM9 molecules. (Schütt et al., 2018)

To showcase the value of our trained KRR model, we apply it to a dataset of 10k diastereomers of parent C7H10O2 isomers. This dataset contains molecular structures, but no HOMO energies. Computing the HOMO energies with DFT would take considerable effort and time. With our KRR model – trained on 32k QM9 molecules represented by the MBTR – we gain an immediate overview of the HOMO energies that occur in the dataset. A histogram of all predicted HOMO energies is shown in Fig. 12. We can see that they are uniformly distributed between -6.8 and -3.2 eV for all diastereomers. The energetic scan allows us to quickly detect molecules of interest in a large collection of compounds. Individual molecules can be easily identified, for instance, the molecule with lowest HOMO energy, molecules with highest HOMO energy or those molecules with average HOMO energies. Various molecules of interest, e.g., structures with HOMO energies in a particular region, could subsequently be further investigated with first-principle methods or experiments to determine their functionality for certain applications. In this fashion, fast energy predictions of our KRR model can be analyzed for structures with desired HOMO energy.

VII Conclusion

In this study, we have trained and tested KRR models on three molecular datasets with different chemical composition to predict molecular HOMO energies. Our comparison between two different molecular descriptors, the CM and the MBTR, shows that the MBTR outperforms the CM on both OE and QM9 due to its ability to encode complex information from the molecular structures. For AA we could find no significant difference in performance between the two representations.

Our work demonstrates that the predictive performance of KRR inherently depends on the complexity of the dataset it is applied to, in addition to the training set size and descriptor. Rapidly decreasing learning curves and low MAEs are achieved for QM9, which is known as a standard benchmark set for ML in molecular chemistry, containing pharmaceutically relevant compounds with rather simple bonding patterns. The same is true for AA which consists of a primitive and restricted collection of amino acids and peptides. The OE dataset, however, comprises large opto-electronic molecules with complicated electronic structures and unconventional functional groups, is much more difficult to learn. It yields almost flat learning curves and considerably higher MAEs. To further improve the predictive power for molecules of technological interest, such as the ones in the OE set, future work should focus on generating larger datasets, devising better descriptors or more sophisticated machine learning methods.

Acknowledgements.

We greatfully acknowledge the CSC-IT Center for Science, Finland, and the Aalto Science-IT project for generous computational resources. This study has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 676580 with The Novel Materials Discovery (NOMAD) Laboratory, a European Center of Excellence and from the Magnus Ehrnrooth Foundation. This work was furthermore supported by the Academy of Finland through its Centres of Excellence Programme 2015-2017 under project number 284621.

Appendix A Chemical names of molecules

Figure 13 shows example molecules from QM9, AA and OE and their chemical names.

Appendix B Choice of kernels and MBTR terms

In Fig. 14, we show learning curves for different parameter configurations of the CM and the MBTR. MAEs correspond to out-of-sample predictions made in a single experiment. We generated learning curves using the MBTR with one-, two- and three-body terms (k=[1,2,3]), two-and three-body terms (k=[2,3]) and only two-body-terms (k=2). For QM9 and OE, we the MBTR with k=[2,3] performs slightly better than the MBTR with k=[1,2,3], while for AA, the performance of k=[2,3] and k=[1,2,3] is equal for larger training set sizes. Therefore, we chose to employ the MBTR with two- and three body terms (k=[2,3]) in this study. Moreover, Fig. 14 reveals that, when the MBTR is used as molecular descriptor, the Gaussian kernel (GK) yields better results than the Laplacian kernel (LK). For the CM, on the other hand, the Laplacian kernel works better than the Gaussian kernel due to discontinuities of the sorted CM.

Appendix C Cross-validated grid search

For MBTR and KRR hyperparameter selection, we employed a 5-fold cross-validated grid search. In a 5-fold cross validation, a given original training set (1k, 2k, 4k, 8k, 16k or 32k) is shuffled randomly and split into 5 equally sized groups. One group (20% of the training set) is taken as a hold-out set for validation. The remaining groups (80% of the training set) are taken as training data. Then, a grid search is performed: Each possible combination of KRR and MBTR hyperparameter values is trained on the training data and evaluated on the validation data.

The assignment of the 5 groups into training data and validation data is repeated 5 times, until each group was used as validation data exactly once and used as training data exactly 4 times. As we repeat the process 5 times, we get 5 MAEs for each possible set of hyperparameters. We consider the average value over these 5 MAEs for each set of hyperparameters and choose the set with lowest average MAE. With the chosen set of optimal hyperparameters, we train the KRR model on the entire original training set (=training data$+$validation data). Finally, we evaluate the trained model on the test set of 10k out-of-sample molecules and report the MAE on the test set in the final learning curve for the original training set size.