# Obtaining transferable chemical insight from solving machine-learning
classification problems: Thermodynamical properties prediction, atomic
composition as good as Coulomb matrix

Leon Alday-Toledo1 Roberto Bernal-Jaquez2 Saul Zapotecas-Martinez3 Jose L.
Mendoza-Cortes4

(1Posgrado en Ciencias Naturales e Ingeniería, Universidad Autónoma
Metropolitana Unidad Cuajimalpa, Ciudad de México, México
2Departamento de Matemáticas Aplicadas y Sistemas, Universidad Autónoma
Metropolitana Unidad Cuajimalpa, Ciudad de México, México
3Coordinación de Ciencias Computacionales, Instituto Nacional de Astrofísica
Óptica y Electrónica, Tonantzintla, Puebla 72840, México
4Department of Chemical Engineering & Materials Science, Michigan State
University, East Lansing, Michigan 48824, United States
<EMAIL_ADDRESS>)

###### Abstract

Machine learning (ML) can be used to construct surrogate models for the fast
prediction of a property of interest. ML can thus be applied to chemical
projects, where the usual experimentation or calculation techniques can take
hours or days for just one sample. In this manner, the most promising
candidate samples could be extracted from an extensive database and subjected
to further in-depth analysis.

Despite their broad applicability, it can be challenging to apply ML methods
to a given chemical problem since a multitude of design decisions must be
made, such as the molecular descriptor to use or the optimizer to train the
model.

Here we present a methodology for the meaningful exploration of a given
molecular problem through classification experiments. This conceptually simple
methodology results in transferable insight on the selected problem and can be
used as a platform from which prediction difficulty is estimated, molecular
representations are tested and refined, and more precise or ambitious projects
can be undertaken. Physicochemical insight can also be obtained.

This methodology is illustrated through the use of multiple molecular
descriptors for the prediction of enthalpy, Gibbs’ free energy, zero-point
vibrational energy, and constant-volume calorific capacity of the molecules
from the public database QM9 [1] with $133\,885$ organic molecules. A
noteworthy result is that for the classification problem we propose, the low-
resolution descriptor ‘atomic composition’ [2] can reach a classification rate
almost on par with the high-resolution ‘sorted Coulomb matrix’ [3, 4, 5]
($>90\%$), provided that an appropriate optimizer is used during training.

###### Contents

  1. 1 Introduction
    1. 1.1 Choice of molecular representation
  2. 2 Theoretical methods
    1. 2.1 An in-depth look: atomic composition
    2. 2.2 An in-depth look: Coulomb-matrix descriptors
    3. 2.3 Choice of the model’s outputs
    4. 2.4 Data preparation
    5. 2.5 Underlying machinery of the neural networks
    6. 2.6 Choice of neural network architecture
  3. 3 Results and discussion
    1. 3.1 Enthalpy
    2. 3.2 Gibbs’ free energy
    3. 3.3 Zero-point vibrational energy
    4. 3.4 Constant-volume heat capacity
  4. 4 Conclusions

##  1 Introduction

Chemical machine learning (ML) seeks to obtain computationally expensive
molecular properties by means of training a predictor with computationally
simpler information. The predictor learns the underlying patterns in the data
and creates a model that connects the inputs with the outputs. If the input
properties are deeply related to the outputs, we expect the predictions to be
accurate.

Chemical ML projects have become feasible in recent years given the broad
availability of modern computational resources, both in terms of hardware and
software. Furthermore, public molecule datasets have emerged, some with
cardinalities in the range of hundreds-of-thousands [1] or even near-billion
and hundreds-of-billions [6, 7].

Molecules have to be encoded by means of a _molecular representation_ or
_descriptor_ in order to be compatible with ML models. This is an underlying
component of any ML model, with a strong influence on the prediction quality.
One project that explores the relationship between descriptor choice and
prediction quality is that of Faber _et al._ [8], who tested nine
representations paired with six regressors for nine target properties. Their
findings are promissory: Several of their models closely follow B3LYP data,
and they could potentially outperform it if a dataset with a higher theory
level were available.

This problem can also be approached from the point of view of mathematics,
treating minimization of prediction error as an optimization problem. One such
example is that of Browning _et al._ [9], who opted to select their training
set in a guided manner, using genetic algorithms. The resulting training sets
can train a NN to a markedly better performance compared to a NN trained with
a randomly-selected training set, which is the usual method. This method was
tested on, and proved true, for a multitude of molecular properties, including
thermodynamical properties such as $H$, $G$, $C_{V}$, electronic properties
such as $E_{\text{HOMO}}$, $E_{\text{LUMO}}$, $E_{\text{gap}}$, and others
such as dipole moment and isotropic polarizability.

One manner to tackle this problem from the chemist’s angle is to consider
other molecular representations that could supply a predictor with new
information and thus improve its performance. However, carrying this endeavor
to completion can prove to be deceptively difficult, as the decision of the
set of candidate descriptors is just one of the many choices that have to be
made before an experiment can be realized. For example, multiple ML models can
be used for prediction, a variety of parameter optimization methods can be
used for training the model, and numerous programming languages and ML
packages can be considered.

In this paper we shall use neural networks (NNs) to model our predictors: they
are a popular choice in terms of literature and implementation, and can be
used for problems such as regression and classification. Their broad
applicability is due to their ‘universal approximator’ character [10, 11].
Indeed, as stated by Hornik _et al._ [11],

> …any lack of success in applications must arise from inadequate learning,
> insufficient numbers of hidden units or the lack of a deterministic
> relationship between input and target.

Note that other ML models or techniques can be applied to chemical problems,
such as support-vector regression ($\varepsilon$-SVR), $\nu$-SVR, and random-
forest regression. These three were tested by Aldosari _et al._ [12] for the
prediction of constant-pressure calorific capacity and entropy of
hydrocarbons. A similar endeavor is that of Liu _et al._ [13], where $PVT$
properties of pure compounds H2O, CO2, H2, and ternary mixtures of them, are
predicted by means of $\nu$-SVR rather than through the development of complex
equations of state. This ML approach yields ‘extremely satisfactory mapping
and prediction results’. Other models such as convolutional neural networks
(CNN) have been previously applied to chemical problems. For example, Krüger
_et al._ [14] use CNNs for the prediction of the reduction potential of
quinones.

Given that we have selected neural networks for this project, additional
choices that must be made are the dataset and optimizer to be used in
training, the molecular encoding (the inputs), target properties (outputs),
the NN architecture (the amount of layers, the amount of neurons in each one,
and the activation functions that connect the layers), and the loss function
(the error metric to minimize). Experiments can be performed once these
settings are chosen.

The findings obtained from these experiments, such as the identification of
high-performance representations, elucidate the selected problem through
resulting transferable knowledge. Insight is gained on the involved
properties, and on possible future experiments that possess increased
precision by selecting a stricter error criteria, or increased ambition by
using a dataset with a greater amount and variety of molecules, such as the
GDB family [15, 16, 6, 7, 17], the biggest of which contains 166 billion
($1.66\text{\times}{10}^{11}$) molecules. Alternatively, other ML methods
could also be tested.

As we shall show, our target properties were enthalpy ($H$), Gibbs’ free
energy ($G$), zero-point vibrational energy (ZPVE), and constant-volume
calorific capacity ($C_{V}$). We trained neural networks (NN) with three
hidden layers for molecule classification according to the value of each of
their properties. Our experiments reveal that the coarse-grained ‘atomic
composition’ representation can be sufficient to classify molecules according
to the value of its $H$, $G$, or ZPVE, whereas $C_{V}$ prediction is
significantly more difficult. Additionally, AdaBelief [18] is shown to be a
noticeably better optimizer choice than Adam [19] for the training of low-
resolution models.

Our dataset of choice is QM9 [1], a popular and public resource constituted of
$133\,885$ organic molecules with up to nine heavy atoms (C, N, O, F). The
molecules are presented as simple text files containing the XYZ coordinates of
the optimized structure, plus additional information such as its harmonic
frequencies, dipole moment norm, $E_{\text{HOMO}}$, $E_{\text{LUMO}}$, and
thermodynamical properties such as zero-point vibrational energy (ZPVE),
constant-volume caloriric capacity $C_{V}$, enthalpy, free energy, and
internal energy. Geometry optimization and property calculations were
performed at the level B3LYP/6-31G(2df,p). We opted for using this dataset
given its popularity and the wealth of properties it provides.

An alternative dataset choice is PC9 [20, 21]. Its molecules also have up to 9
heavy atoms, but they show a greater diversity regarding bond distances and
functional groups. PC9 is constituted by $99\,234$ molecules, $80\,877$ of
which are not present in QM9; these include $4442$ radicals and $883$
triplets; compare this with QM9: all its molecules are closed-shell and
neutral. The entries in the PC9 database contain the molecular geometries,
atomic charges, total energies, and those of $\text{HOMO}-1$, HOMO, LUMO,
$\text{LUMO}+1$, and were calculated at the theory level B3LYP/6-31G(d).

We used $k$-fold cross-validation with $k=5$, therefore all error decay plots
and classification tables show the average values for five experiments with
different training and test sets.

Additionally, we modelled our neural networks using the scientific-computing-
oriented julia programming language [22], and the elegant and extensible
machine-learning package Flux.jl [23, 24]; other julia packages were also
used[25, 26, 27, 28, 29, 30, 31, 32, 33].

###  1.1 Choice of molecular representation

In this paper, we shall test multiple molecular representations and compare
and analyze their performance so as to identify their relative performances
and possible beneficial augmentations. Three representations we term
_standalone_ , since these use a sole family of properties as NN input. Two
additional representations are _hybrid_ or _mixed_ , since they are the
concatenation of two representations:

  * •

Two standalone representations are highly traditional: the sorted Coulomb
matrix (SCM) and its eigenvalues (EV) [3, 4, 5].

  * •

The third one is the atomic composition (AC) [2].

  * •

The two hybrids are the concatenations of SCM and AC, and of EV and AC. We
abbreviate them SCMAC and EVAC, respectively.

The SCM and SCMAC representations encode the molecular stoichiometry and
geometry, and their sizes grow quadratically with respect to the molecule’s
atom count. Therefore, we consider these _high-resolution_ (or fine-grained)
representations. On the other hand, the EV and EVAC representations scale
linearly with the system’s size, and AC has constant size for any particular
dataset, regardless of the selected molecule. Therefore, these three
representations are considered to be _low-resolution_ (or coarse-grained).

##  2 Theoretical methods

###  2.1 An in-depth look: atomic composition

The AC representation lists the amount of atoms in a molecule, considering not
the molecule alone but the atom types featured across the entire dataset:
Formaldehyde, H2CO, is QM9’s sixth molecule. Since QM9 includes the atom types
$(\ch{H},\ch{C},\ch{N},\ch{O},\ch{F})$, formaldehyde’s AC representation
within the QM9 context is the vector $(2,1,0,1,0)$.

It can be seen that AC is a low-resolution descriptor and it doesn’t represent
molecules uniquely, namely isomers. Nevertheless, we consider it to encode
chemically important information and to be a representation worth testing for
these reasons:

  * •

Burden and Winkler [34, 35] have employed an idea similar to AC for prediction
of molar refractivity, hydrophobicity and biological activity of drug-like
molecules.

  * •

Tchagang and Valdés [2] show that the EV and SCM representations can be
augmented by concatenation of the AC representation. They predict atomization
energy using Bayesian-regularized NNs, for the molecules in the QM7 database
[36]. The augmented models reach markedly lower prediction errors: their
lowest mean absolute error (MAE) is
$3.0\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$, reached by an SCMAC
model, whereas the lowest MAE reached by SCM models is
$4.4\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$.

Their findings show that augmented or mixed representations are a simple yet
effective manner to improve prediction quality. Therefore, AC is a promising
auxiliary representation.

  * •

a wealth of literature exists regarding group-contribution methods. These seek
to approximately calculate certain molecular properties from the structure
alone. The first of its kind was proposed by Riedel [37] in 1949. His method
calculates the critical pressure of organic compounds, based only on two
variables: the molecule’s mass and an additive quantity $\varphi$ calculated
from its structure.

Other group-contribution methods calculate the critical properties of pure
compounds [38, 39, 40, 41, 42, 43], boiling and freezing temperatures [41, 42,
44], constant-pressure calorific capacity [45, 41, 12], entropy [45, 12],
enthalpy and Gibbs’ free energy of vaporization, fusion, and formation [45,
41, 42], and liquid viscosity [41].

There’s also Fredenslund _et al._ ’s UNIFAC method [46], which approximates
activity coefficients for components in binary and ternary nonideal liquid
mixtures.

While many of these methods consider detailed molecular descriptions involving
the contribution of atoms, bonds, and functional groups, we consider that AC
can be tested for exploring the thermodynamical landscape of the QM9 dataset.
AC can be considered to be analogous to Benson’s law of additivity of atomic
properties [45], which is the zero-order approximation to the law of
additivity of molecular properties.

Therefore, one of the objectives of this work is to evaluate the performance
of AC and compare it with the established representations EV and SCM.
Furthermore, we shall also test the ability of mixed representations involving
AC and the Coulomb-matrix-derived representations. Since we study the relation
between a molecule’s stoichiometry and some of its thermodynamical properties,
this work can be considered an ML-informed successor to the aforementioned
group-contribution methods.

###  2.2 An in-depth look: Coulomb-matrix descriptors

The sorted Coulomb matrices were calculated using the molecular geometries and
the python [47] package QML [48]. The eigenvalues of the matrices were then
calculated with the package numpy [49].

The Coulomb matrix (CM) representation [3] is a conceptually simple encoding
of the geometry. Its elements are defined

$C_{ij}=\left\\{\begin{array}[]{cc}\dfrac{1}{2}Z_{i}^{2.4}&\text{ if
}i=j\\\\[8.00003pt]
\dfrac{Z_{i}Z_{j}}{\lVert\vec{R}_{i}-\vec{R}_{j}\rVert}&\text{ if }i\neq j\\\
\end{array}\right.$ (1)

where $Z_{i}$ is the nuclear charge of atom $i$, and $\vec{R}_{i}$ refers to
its Cartesian coordinates in bohr [50]. The formula for the diagonal elements
was obtained from fitting the potential energies of the free atoms to their
nuclear charges; the non-diagonal elements represent the Coulomb repulsion
between every pair of atoms [3].

The CM representation is invariant to translation or rotation, but not to atom
reordering in the coordinate list. On the other hand, the eigenvalue
representation is invariant to these three operations, and, despite its low
dimensionality, Rupp _et al._ [3] show it can be used to predict atomization
energies, outperforming other methods such as bond counting [51] or
semiempirical PM6 [52].

Furthermore, Coulomb matrices can be made invariant to atom reordering by
means of sorting its rows and columns [4, 5]: the rows and columns of an
unsorted CM may be uniquely sorted by arranging them according to sorting the
nondecreasing order of their norms. Thus, the sorted Coulomb matrix
representation (SCM) is obtained, which we use in this paper. Since the SCM is
a symmetrical matrix, we only need to consider its nonrepeated elements for
training.

###  2.3 Choice of the model’s outputs

The choice of a model’s output is deeply related to the input features and the
model’s purpose: if the selected inputs are numerous, then an adequate model
would need a high amount of hidden neurons and layers. This is exacerbated if
the model is tasked with identifying or approximating multiple sample
properties, or if the relation between the inputs and outputs is complex.

In the interest of keeping the model’s training affordable, and our results
simple enough to be analyzed, we shall perform classification experiments,
rather than regression ones. In the latter, the model outputs a continuous
variable signifying the predicted value of a given property of a sample. A
target test-set average prediction error can then be defined, and the models’
performance can be compared to it. On the other hand, in a classification
experiment, the metric to optimize is the classification rate.

Therefore, we must define the classes for our experiment: we can define non-
overlapping intervals such that a molecule is considered to belong to class
$i$ if its value of a certain property is within the $i$-th interval. These
intervals needn’t be equally-sized, they merely need to be significantly
populated, since having great disparity in the cardinality of the classes can
lead to heterogeneous prediction behavior and estimating the generalization
error can become difficult.

The advantages of this classification technique are that both the model’s
global and per-class performance are evaluated. The latter can reveal the
worst-performing regions of the dataset; focusing on them to improve the
model’s overall performance can be more effective than focusing on improving
the global behavior. Furthermore, once our classification model is trained, it
can be used to swiftly filter a database and find the subset of molecules
whose value of a given property is within a desired range. Afterwards, these
top candidates can be subjected to computationally expensive electronic
structure calculations in order to narrow down the list of candidate
structures even further.

On the other hand, if our experiments were unable to find a well-performing
model, some of the obtained knowledge can be transferable to future
experiments, such as the classes that show the easiest prediction, which
representations among a selection have the best performance and how they could
be improved, or which optimizer is the best suited for training models with
the selected inputs and outputs. Therefore, this exploration would still
inform and support any other future experiments.

Some classification experiments could have classes with widths beyond the
typical computational-chemistry benchmarks. Note, however, that as long as the
cardinalities of the classes are not significantly heterogeneous, there is no
formal restriction in terms of how many or how thin the classes of a given
experiment may be. Therefore, subsequent experiments could use greater
granularity by increasing the amount of classes (and thus decreasing their
width); regression experiments would also be an option.

If we were to seek a regression model that predicts the value of one property,
the NN’s output layer would have a single neuron containing the predicted
value for a particular sample. On the other hand, for our classification
experiments we have as many output neurons as there are classes: neuron $i$
contains the degree of similarity between the current sample and class $i$.
Therefore, by deciding how many classes or intervals of interest our
experiment shall have, we are also deciding the amount of output neurons in
the model.

###  2.4 Data preparation

We standardized the values of the inputs and output properties before
performing the experiments, as is commonly done in ML experiments so that
inputs and outputs with different absolute values can be used in the same
model.

The molecular descriptors we selected have shown to be effective for the
prediction of thermodynamical properties. Therefore, the output properties we
selected are enthalpy ($H$), Gibbs’ free energy ($G$), zero-point vibrational
energy (ZPVE), and constant-volume calorific capacity ($C_{V}$). Some of their
distribution properties are shown in table 1.

We must define the amount and thresholds of the classes for each one of these
properties. For our experiments, we shall divide the standardized outputs in
25 bins, to be grouped into classes with comparable cardinalities. This leads
to $H$ and $G$ having five classes; ZPVE and $C_{V}$ having seven; subsequent
experiments could use more bins, this leading to more and thinner classes. The
current bins and classes are shown in figures 1 to 4. See also tables 2 and 3
for their thresholds and populations.

These were the lengths of the input representations we used:

  * •

AC: as mentioned above, this representation is of length $5$.

  * •

EV: since the largest QM9 molecules are constituted by $29$ atoms (hydrogens
included), our EV representation was fixed at length $29$.

  * •

SCM: our matrices are of size $29\times 29$, with a total of $29^{2}=841$
properties. However, the duplicate nondiagonal elements were discarded, and
thus this representation had length $\frac{29\cdot 30}{2}=435$. For all
molecules with less than 29 atoms, their SCMs were padded with rows and
columns of zeroes until size $29\times 29$ was reached, as suggested by the
method’s authors [3, 4].

  * •

EVAC and SCMAC: These representations have respective lengths $34$ and $440$.
These keep EV and SCM’s ability to represent molecules uniquely, unlike AC.

| min | max | $\bar{x}$ | $\sigma$
---|---|---|---|---
$H$/$\mathrm{\text{$E$}_{\mathrm{h}}}$ | $-714.56$ | $-40.48$ | $-411.53$ | $40.06$
$G$/$\mathrm{\text{$E$}_{\mathrm{h}}}$ | $-714.60$ | $-40.49$ | $-411.58$ | $40.06$
ZPVE/$\mathrm{\text{$E$}_{\mathrm{h}}}$ | $0.0160$ | $0.274$ | $0.149$ | $0.0333$
$C_{V}$/$\mathrm{cal}\text{\,}{\mathrm{mol}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ | $6.002$ | $46.969$ | $31.601$ | $4.062$
Table 1: Statistical parameters of $H$, $G$, ZPVE and $C_{V}$ in the QM9 database [1]. Figure 1: Histogram of $H$ and its designated classes. Figure 2: Histogram of $G$ and its designated classes. Figure 3: Histogram of ZPVE and its designated classes. Figure 4: Histogram of $C_{V}$ and its designated classes. | $H$ | $G$
---|---|---
Class | bins | values/$\mathrm{\text{$E$}_{\mathrm{h}}}$ | proportion | bins | values/$\mathrm{\text{$E$}_{\mathrm{h}}}$ | proportion
$c_{1}$ | $<11$ | $[$-714.559\,209$,$-444.925\,572\,200\,000\,03$)$ | $0.173\,358$ | $<11$ | $[$-714.602\,138$,$-444.960\,721\,6$)$ | $0.173\,358$
$c_{2}$ | $11$ | $[$-444.925\,572\,200\,000\,03$,$-417.962\,208\,52$)$ | $0.319\,281$ | $11$ | $[$-444.960\,721\,6$,$-417.996\,579\,96$)$ | $0.321\,238$
$c_{3}$ | $12$ | $[$-417.962\,208\,52$,$-390.998\,844\,839\,999\,95$)$ | $0.197\,289$ | $12$ | $[$-417.996\,579\,96$,$-391.032\,438\,320\,000\,04$)$ | $0.195\,332$
$c_{4}$ | $13$ | $[$-390.998\,844\,839\,999\,95$,$-364.035\,481\,159\,999\,96$)$ | $0.196\,624$ | $13$ | $[$-391.032\,438\,320\,000\,04$,$-364.068\,296\,68$)$ | $0.196\,624$
$c_{5}$ | $>13$ | $[$-364.035\,481\,159\,999\,96$,$-40.475\,116\,999\,999\,955$)$ | $0.113\,448$ | $>13$ | $[$-364.068\,296\,68$,$-40.498\,597\,000\,000\,02$)$ | $0.113\,441$
Table 2: Thresholds and proportions of the classes for $H$ and $G$. | ZPVE | $C_{V}$
---|---|---
Class | bins | values/$\mathrm{\text{$E$}_{\mathrm{h}}}$ | proportion | bins | values/$\mathrm{cal}\text{\,}{\mathrm{mol}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ | proportion
$c_{1}$ | $<11$ | $[$0.015\,950\,999\,999\,999\,993$,$0.119\,148\,200\,000\,000\,01$)$ | $0.183\,949$ | $<14$ | $[$6.001\,999\,999\,999\,999$,$27.304\,84$)$ | $0.135\,071$
$c_{2}$ | $11$ | $[$0.119\,148\,200\,000\,000\,01$,$0.129\,467\,920\,000\,000\,01$)$ | $0.102\,103$ | $14$ | $[$27.304\,84$,$28.943\,52$)$ | $0.115\,084$
$c_{3}$ | $12$ | $[$0.129\,467\,920\,000\,000\,01$,$0.139\,787\,640\,000\,000\,02$)$ | $0.124\,794$ | $15$ | $[$28.943\,52$,$30.5822$)$ | $0.151\,727$
$c_{4}$ | $13$ | $[$0.139\,787\,640\,000\,000\,02$,$0.150\,107\,360\,000\,000\,02$)$ | $0.122\,725$ | $16$ | $[$30.5822$,$32.220\,88$)$ | $0.164\,567$
$c_{5}$ | $14$ | $[$0.150\,107\,360\,000\,000\,02$,$0.160\,427\,08$)$ | $0.127\,901$ | $17$ | $[$32.220\,88$,$33.859\,56$)$ | $0.150\,345$
$c_{6}$ | $15,16$ | $[$0.160\,427\,08$,$0.181\,066\,52$)$ | $0.165\,276$ | $18$ | $[$33.859\,56$,$35.498\,24$)$ | $0.117\,018$
$c_{7}$ | $>16$ | $[$0.181\,066\,52$,$0.273\,944\,000\,000\,000\,1$]$ | $0.173\,253$ | $>18$ | $[$35.498\,24$,$46.969$]$ | $0.166\,18$
Table 3: Thresholds and proportions of the classes for ZPVE and $C_{V}$.
Figure 5: Deep neural network (DNN) with three hidden layers.

###  2.5 Underlying machinery of the neural networks

An illustrative drawing of our NNs can be found in figure 5. They predict the
properties of a sample through a process known as _forward-propagation_ [53],
in which the molecule’s properties are propagated through the NN’s neurons by
combination of the inputs with the NN’s parameters (its weights and biases),
which are processed by the neural network’s _activation function_ before they
reach the next layer.

For classification experiments, the _sigmoid function_ (eq. 2) is a common
choice. Let $z$ be the model’s bias term plus the sum of the product of the
model’s weights and the sample’s properties. Therefore, when $z\to-\infty$,
$f(z)\to 0$, and when $z\to\infty$, $f(z)\to 1$.

$f(z)=\frac{1}{1+\exp(-z)}$ (2)

At the beginning of an ML experiment, the model’s parameters are initialized
randomly. Since our NNs use the sigmoid activation function, initialization is
done through the Xavier initialization scheme111Originally proposed by Xavier
Glorot and Joshua Bengio [54], this method is also known as Glorot
initialization. [54]. The output layer’s neurons show the similarity between a
sample and each of the classes: a perfect model would give a score $1.0$ for
the classes the sample belongs to; $0.0$ for those it does not.

Training of the ML model is therefore intended to bring the model’s
performance as close as possible to the ideal model’s: the model’s parameters
are iteratively adjusted such that the experiment’s _loss function_ is
minimized, which quantifies the model’s prediction error for all samples in
the training set.

As is typically done in classification experiments, we use the binary cross-
entropy loss function:

$J_{m}(W,b)=-\frac{1}{N}\sum_{i=1}^{N}\left[y_{i}\cdot\ln(a_{i})+(1-y_{i})\ln(1-a_{i})\right]$
(3)

where $J_{m}$ is the error for sample $m$, $N$ is the amount of cells in the
output layer (also the amount of classes in the experiment), $y_{i}$ is the
label of sample $m$ for class $i$, and $a_{i}$ is the model’s prediction for
class $i$ of sample $m$, a continuous value $\in(0,1)$.

Furthermore, the loss minimization process is performed by means of an
optimizer. Many popular choices are gradient-based, such as Adam [19, 55], a
modification of the traditional gradient descent (GD) method. It is an
adaptive method since it has a separate step size (_learning rate_) for each
parameter, which are reduced at different rates depending on their respective
slopes; compare this with GD: it has a fixed, common step size for all
parameters even if some parameters have converged and some others have not.

Adam builds upon two gradient-based optimization methods: AdaGrad [56] and
RMSProp [57]. It possesses the former’s ability to deal with sparse gradients
and the latter’s ability to work with non-stationary problems. We use it for
our experiments given its popularity and performance.

We selected a second optimizer: AdaBelief [18], a recent variant of Adam. It
seeks to maintain the adaptive methods’ stable and fast training, while
attaining a test-set (generalization) error similar to that of SGD, which can
outperform adaptive methods [58]. One of its distinct components is the
_belief_ that “the exponential moving average (EMA) of the noisy gradient is a
prediction of the gradient at the next time step”. Therefore, it compares the
current EMA with the observed gradient for the next step. If they differ, the
observed gradient is distrusted and parameters are updated with a smaller step
size. Conversely, if they are similar, the effective step size is larger.
Given its recency and solid theoretical foundations, we decided to include it
our project.

We used the default Adam and Adabelief hyperparameters, just as shown by the
authors in their original papers: $\eta=0.001$, $\beta_{1}=0.9$,
$\beta_{2}=0.999$, $\epsilon=$1\text{\times}{10}^{-8}$$.

###  2.6 Choice of neural network architecture

For any given problem, the choice of the amount of layers and their amount of
nodes (the model’s architecture) is a nontrivial question that requires manual
experimentation and testing of multiple candidate architectures; while no
certain formulas or procedures exist for the analytical calculation of a
perfect architecture, this search can be eased if domain knowledge is
available. Another approach to ease the architecture search is through the
discipline of ‘neural network architecture search’ (NAS), which encompasses
multiple methods that seek to automate this search [59, 60].

Some useful guidelines exist, which can be used to obtain candidate
architectures intuitively: a neural network should have a funnel-like shape,
with early hidden layers having higher neuron counts than the latter layers,
and no sharp cardinality drops between one layer and another. Also, the
prediction and error performances should be similar for the training and test
sets; this indicates the current architecture is experiencing no significant
underfit or overfit, and thus can be considered an adequate architecture for
the task at hand.

Our input vectors have lengths from $5$ to $440$. Such a wide range warrants
using multiple architectures, depending on the amount of features, rather than
a single common architecture. Given the aforementioned non-triviality of the
neural architecture search, we tested multiple neural network architectures
until settling for three-hidden-layer models with the following cardinalities:

  * •

For AC, inp:12:10:8:out,

  * •

for EV and EVAC, inp:26:20:14:out,

  * •

for SCM and SCMAC, inp:240:120:60:out.

where inp is the aforementioned length of the respective feature vector (the
standalone representation’s lengths are $5$, $29$, and $435$), and out is the
amount of classes of the target property ($5$ for $H$ and $G$, $7$ for ZPVE
and $C_{V}$).

##  3 Results and discussion

###  3.1 Enthalpy

Figure 6 shows the decay of the training-set error (binary cross-entropy
error, equation 3) after each epoch of training. While a low training error is
desirable, a model’s ultimate test is the _generalization_ or _out-of-sample_
error: how would our model perform for heretofore unknown samples, such as
those from a future dataset? A simple way to approximate this value is by
measuring the model successful classification rate for the molecules in the
test set, which didn’t take part in the training. This is shown in table 4. We
consider a sample to be classified correctly if the neuron of the sample’s
class has the biggest score among the output neurons.

As can be seen in figure 6, the three AdaBelief curves reach near-zero values,
whereas the Adam curves don’t. This is consistent with table 4, where
AdaBelief is largely unconcerned by the input features; Adam can reach $>99\%$
classification performances for the high-resolution inputs, while the low-
resolution scores are $\in(0.57,0.88)$. Therefore, the optimizer choice is
especially important for low-resolution models: the correct optimizer choice
can allow a low-resolution model to reach nearly indistinguishable end-point
training error and test-set classification performance as a high-resolution
model. Furthermore, note that while reaching this performance level requires a
training with a greater amount of epochs than that of the high-resolution
models, these epochs come at a lower computational cost, given the stark
difference in the models’ overall amount of neurons.

AC is shown to be a capable standalone descriptor for $H$ prediction. Its
influence in the performance of the mixed-representation models is also
positive: figure 6 shows how the mixed-representation models generally train
faster than the standalone EV and SCM ones. While the final error of the four
models is practically identical, training could be stopped after the cost
function has converged, thus the mixed-representation model would spend less
computer time than the standalone.

We can draw further noteworthy observations by considering the least-
performant experiments:

  * •

The EVAC representation is the richest out of the low-resolution inputs, given
that it includes both EV and AC. Therefore, its performance is the best among
its peers.

  * •

Consider table 4. Even if we didn’t have access to the AdaBelief optimizer and
its consistently high performances, we could still perform training
experiments with Adam, and improve their performances through the use of
richer descriptors.

  * •

regarding the Adam models, the classification performance of EV and AC are,
respectively, $0.5754$ and $0.7279$, compared to EVAC’s $0.8783$. We can thus
make an argument similar to Hammond’s postulate [61], which is usually
expressed as the following corollary: _for a given transition state (TS) that
connects two minima structures, the TS exhibits greater geometrical similarity
to the minimum its energy is closest to_. In our ML context, we can restate it
in the following manner: _let A be a standalone representation and M a mixed
representation that includes A. Therefore, if A’s prediction performance is
closest, among its peers, to M’s, then A is the most influential
representation in M’s, among all the representations that constitute M_. This
is to say, A is the standalone representation that provides the most valuable
information among those in M.

In terms of enthalpy prediction according to our classification experiment, we
can thus conclude that AC is a bigger contributor to EVAC’s high prediction
quality than EV.

This is a simple procedure that can provide us with valuable insight, and we
term it _Hammond-like analysis_.

Figure 6: Training-set error decay of the enthalpy models. Input |  |
---|---|---
properties | Adam | AdaBelief
AC | 0.7279 | 0.9921
EV | 0.5754 | 0.9915
EVAC | 0.8783 | 0.9939
SCM | 0.9943 | 0.9947
SCMAC | 0.9948 | 0.9953
Table 4: Global classification performance of the enthalpy models for test-set
molecules.

###  3.2 Gibbs’ free energy

This property shows a very similar behavior to enthalpy. Once again, the high-
resolution inputs achieve near-perfect performance; the low-resolution models
have an excellent behavior when trained with AdaBelief ($>99\%$ classification
performance, table 5), whereas the Adam-trained models all have scores
$\in(0.75,0.94)$.

One noteworthy difference between the $H$ and $G$ models is that while the AC-
Adam performance remains steady at $73\%\sim 76\%$ for both properties, EV-
Adam goes from $0.5754$ for $H$ prediction to $0.8987$ for $G$. This suggests
the eigenvalues are more closely related to $G$ than they are to $H$. The EVAC
representation is therefore also better suited for $G$ prediction than to $H$;
this is visible in tables 4 and 5: the Adam scores are, respectively, $0.8783$
and $0.9359$. AdaBelief exhibits minute shifts when comparing $H$ to $G$, but
they are consistent with the aforementioned trend.

While the _a priori_ table 2 and distribution plots 1 and 2 show that $H$ and
$G$ are nearly identical, the _a posteriori_ training curve plots 6 and 7 and
table 4 and 5 show that these properties are distinct; the Hammond-like
analysis shows that EV is more important for $G$ prediction than AC, whereas
the converse was true for $H$.

Figure 7: Training-set error decay of the Gibbs’ free energy models. Input |  |
---|---|---
properties | Adam | AdaBelief
AC | 0.7580 | 0.9913
EV | 0.8987 | 0.9932
EVAC | 0.9359 | 0.9952
SCM | 0.9954 | 0.9955
SCMAC | 0.9959 | 0.9959
Table 5: Global classification performance of the Gibbs’ free energy models
for test-set molecules.

###  3.3 Zero-point vibrational energy

This property has a qualitatively similar behavior to that of $H$ and $G$.
Once again, the high-resolution descriptors have the best performances, while
the low-resolution descriptors perform far better with AdaBelief than with
Adam. Furthermore, the Hammond-like analysis suggests that the AC
representation is a bigger contributor to the EVAC performance than EV, and it
also contains more readily-available relevant information for ZPVE prediction.

Figure 8 shows this property’s increased difficulty: the low-resolution models
no longer reach near-zero error values; they only reach $<0.1$ on two out of
six decays. High-resolution descriptions’ performance is also worsened,
reaching $0.0221$ compared to $H$’s $0.0041$ and $G$’s $0.0033$.

Let us discuss the differences between ZPVE and the two previous properties:

  * •

AdaBelief reached classification scores $>99\%$ for all $H$ and $G$ models,
whereas the ZPVE scores are $\in(0.43,0.97)$. If the EV result is ignored,
then the scores are $\in(0.91,0.97)$, which is nonetheless a greater
dispersion than that of $H$ and $G$.

  * •

In order to better understand this property’s comparatively high difficulty,
table 6 includes two additional columns, showing AdaBelief’s classification
performance for the two classes that exhibit the lowest overall scores.

This per-class analysis shows where our model-refinement efforts may be best
employed: rather than seeking a different descriptor with a better global
score, we could instead search for a representation that performs adequately
for $c_{5}$ and $c_{6}$. These additional features could be concatenated to
any of the five aforementioned molecular descriptions, and we expect this new
mixed representation to show higher $c_{5}$, $c_{6}$, and global scores.

EV can be considered such a representation to AC: its $c_{6}$ performance is
better than its dataset average, and its inclusion allows EVAC to nearly rival
SCM and SCMAC’s performance, both for $c_{6}$ and globally. On the other hand,
EV’s contribution to EVAC’s $c_{5}$ performance is much less significant.

Figure 8: Training-set error decay of the ZPVE models. Input | Global | AdaBelief
---|---|---
properties | Adam | AdaBelief | $c_{5}$ | $c_{6}$
AC | 0.5069 | 0.9106 | 0.8878 | 0.6686
EV | 0.3460 | 0.4302 | 0.2830 | 0.4717
EVAC | 0.7867 | 0.9480 | 0.8959 | 0.8821
SCM | 0.9573 | 0.9550 | 0.9161 | 0.9070
SCMAC | 0.9631 | 0.9644 | 0.9307 | 0.9349
Table 6: Classification performance of zero-point vibrational energy for test-
set molecules.

###  3.4 Constant-volume heat capacity

Let us first focus on the similarities between this property and the previous
three: AdaBelief outperforms Adam for all models, especially in the low-
resolution experiments. Furthermore, the Hammond-like analysis suggests that
AC is a better descriptor for $C_{V}$ than EV.

In terms of differences, however, $C_{V}$ shows the hardest prediction
difficulty. As shown in figure 9, none of the low-resolution models reach
errors below $0.2$, whereas the high-resolution models remain $>0.15$ within
2000 epochs. Consequently, the classification scores (table 7) are
$\in(0.25,0.74)$.

Another difference is visible under AdaBelief: EVAC’s classification score is
well behind that of the high-resolution models; the previous differences were
all $<1\%$.

In order to shed some light on the unique behavior of $C_{V}$, let us turn to
table 8. This table shows the global and per-class classification performance
of AdaBelief alone. This is a more granular display than that of global-only
tables, and it provides us with greater insight on the unique behavior of this
property.

We can see that AC and EV reach their highest scores for the extrema classes
$c_{1}$ and $c_{7}$ ($>68\%$), whereas the intermediate classes have scores
$<0.5$, some even reaching $<0.1$. That said, the EVAC model shows an improved
behavior across all classes; the extrema scores almost reach SCM and SCMAC’s,
most intermediate classes become $>0.51$, while the exception $c_{6}$ still
shows a marked improvement at $0.4398$. This shows that the molecules of the
extrema classes have $C_{V}$ values more closely linked to their atomic
composition or eigenvalue representation, whereas the molecules of the
intermediate classes have other factors influence their $C_{V}$.

The low-resolution scores for $c_{2}$ are also noteworthy: while both
standalones have $<0.1$ scores, their mixed representation has performance on
par with the other intermediate classes; this sudden performance increase is
one of the major reasons behind EVAC’s improved global score. Therefore,
extrapolating the behavior of mixed representations by comparison of the
standalone representations alone can lead to synergistic sets of input
features remaining undiscovered. $c_{6}$ shows a similar behavior.

The performance of any one of the aforementioned representations, low-
resolution especially, can be improved significantly by addition of another
descriptor with a high score for the intermediate classes. This is especially
true for $c_{6}$, where the gap between low- and high-resolution descriptors
is the widest. By improving this class’ performance, the cheaper
representations would once again approach the more expensive ones.

Alternatively, any of the five representations above can be used to enrich
another representation, if its extrema classes’ performance is to be improved.
In general, a similar procedure can be performed for several candidate
representations: once their individual capabilities are found, appropriate
high-performance mixed representations can be created. Physical or chemical
insight would also be obtained.

Figure 9: Training-set error decay of the $C_{V}$ models. Input |  |
---|---|---
properties | Adam | AdaBelief
AC | 0.3217 | 0.4174
EV | 0.2548 | 0.3683
EVAC | 0.5287 | 0.6185
SCM | 0.7143 | 0.7213
SCMAC | 0.7285 | 0.7371
Table 7: Global classification performance of $C_{V}$ for test-set molecules. Input | AdaBelief
---|---
properties | global | $c_{1}$ | $c_{2}$ | $c_{3}$ | $c_{4}$ | $c_{5}$ | $c_{6}$ | $c_{7}$
AC | 0.4174 | 0.6884 | 0.0879 | 0.4956 | 0.3786 | 0.3423 | 0.0178 | 0.7424
EV | 0.3683 | 0.7422 | 0.0367 | 0.2732 | 0.2274 | 0.3816 | 0.0000 | 0.7680
EVAC | 0.6185 | 0.8349 | 0.5670 | 0.5500 | 0.5363 | 0.5156 | 0.4398 | 0.8409
SCM | 0.7213 | 0.8734 | 0.6768 | 0.6676 | 0.6511 | 0.6489 | 0.6320 | 0.8755
SCMAC | 0.7371 | 0.8804 | 0.6856 | 0.6871 | 0.6828 | 0.6598 | 0.6374 | 0.8963
Table 8: Per-class performance of the constant-volume heat capacity models
trained with AdaBelief, for classification of test set molecules.

##  4 Conclusions

We have tested the performance of 3-hidden-layer NNs for the classification of
organic molecules with up to nine heavy atoms, according to the value of their
$H$, $G$, ZPVE, and $C_{V}$, using coarse-grained experiments.

When low-resolution molecular descriptors are used to train a NN, $>99\%$
classification performance can be reached for $H$ and $G$, both with five
classes. ZPVE was divided in 7 classes, and can be classified with $>90\%$
rates. Note that the outcome depends on the optimizer: AdaBelief routinely
achieves high classification performance for most descriptors; the popular
Adam optimizer is unable to reach such performance for low-resolution inputs.

These results show that, depending on the precision target, low-resolution
descriptors with smaller architectures and appropriate optimizers can be
sufficient; expensive experiments with high-resolution features and a high
neuron- and layer-count may be forgone entirely.

Alternatively, low-resolution chemically-significant representations such as
AC may contain highly influential descriptors with new readily-available
information that can be extracted from the inputs directly, rather than
learned or derived by the NN during training. Therefore, aggregating these
descriptors into a mixed representation results in a more descriptive
molecular representation, which can lead to improved model performance and
reduced training time. This holds true even if no ideal optimizer is known,
and can allow one to undertake future experiments seeking more ambitious and
chemically diverse datasets or tighter precision targets.

High-resolution representations were also tested. For these experiments, the
optimizer choice is much less influential in the outcome: Adam and AdaBelief
have very similar performances.

Despite their higher training cost, the performance of the high-resolution
models is only slightly ahead of the low-resolution ones, for all but one
output property: $C_{V}$, which shows the highest prediction difficulty for
all optimizers and input descriptors; the best-performant $C_{V}$ model uses
AdaBelief and the SCMAC representation, reaching a test-set classification
score of $0.7371$. Thus, $C_{V}$ prediction can still be improved.

As we have shown previously, one way to improve the prediction quality for the
more difficult properties is through the creation of an appropriate mixed
representation. The per-class and Hammond-like analyses, while simple, can
nevertheless provide very valuable information and guide the construction of
such a representation. The former can be more revealing than a global-
performance analysis in identifying the current descriptor’s underperforming
classes; the latter can reveal the relative influence of multiple descriptors
that take part in a given mixed representation. For example, the Hammond-like
analysis on $H$ and $G$ shows that the eigenvalue representation is more
significant for the prediction of $G$.

While the Adam optimizer is typically behind AdaBelief, it remains a valuable
tool since it can help us understand whether a given molecular descriptor
presents more learning-ready information for the prediction of a given output
property. Such findings can be used for engineering larger-scale experiments,
with more descriptors and molecules.

The classification approach presented in this paper can provide information of
the mathematical or physicochemical kind regarding the selected training
optimizer, inputs, outputs, and the relations between one another. This
insight can be used to guide and hone further experiments showing more
ambition. For example, a different dataset could be used, such as PC9 [20, 21]
given its greater chemical diversity while remaining of similar size to QM9
[1]; the GDB family [15, 16, 6, 7] is another interesting option, which lists
molecules with up to 17 heavy atoms.

Alternatively, experiments with a higher precision goal could be undertaken,
perhaps now using narrower intervals of interest, or using a greater amount of
bins, thus resulting in a greater amount of classes while maintaining the
relative similarity between the cardinality of each class. This applies to
$H$, $G$, and ZPVE, given their positive performance for our experiments with
25 bins and 5 or 7 classes. Regression experiments could also be explored.

For $C_{V}$, on the other hand, prediction quality must be refined before more
granular experiments are to be attempted. Let us consider the low- and high-
resolution descriptors we used: the main difference between them is the
explicit encoding of interatomic distances, and its inclusion had a positive
effect in the prediction quality. This is consistent with Einstein’s and
Debye’s statistical-mechanics models for $C_{V}$ of solids [62, 63], where
bond strength is an important factor.

Therefore, desirable representations for future testing would include bond
information. Examples of such descriptors are Bag of Bonds [64], bonds-and-
angles machine learning (BAML) [65], histograms of distances, angles, and
dihedral angles [8], Behler-Parrinello or Smith symmetry functions [66, 67,
68]. With the performance baseline established in the aforementioned
experiments, comparisons involving new descriptors are attainable.

## Acknowledgements

The authors thank Dr. Luis Alarcón-Ramos for his very kind help on hardware-
related matters, as well as prompt technical assistance.

The authors thank Valentin Bogad for very swift and in-depth assistance on
diverse matters pertaining implementation and regarding the julia programming
language.

L. A.-T. acknowledges support from scholarship 2020-000013-01NACF-11160 from
CONACyT.

L. A.-T. thanks Dr. M. Davar for persistent discussion and optimization on
diverse writing-adjacent aspects.

## References

  * [1] Raghunathan Ramakrishnan, Pavlo O. Dral, Matthias Rupp and O. Anatole von Lilienfeld  “Quantum chemistry structures and properties of 134 kilo molecules”  In _Scientific Data_ 1, 2014  DOI: 10.1038/sdata.2014.22
  * [2] Alain B. Tchagang and Julio J. Valdés  “Prediction of the Atomization Energy of Molecules Using Coulomb Matrix and Atomic Composition in a Bayesian Regularized Neural Networks”  In _Artificial Neural Networks and Machine Learning – ICANN 2019: Workshop and Special Sessions_ Cham: Springer International Publishing, 2019, pp. 793–803  DOI: 10.1007/978-3-030-30493-5_75
  * [3] Matthias Rupp, Alexandre Tkatchenko, Klaus-Robert Müller and O. Anatole von Lilienfeld  “Fast and accurate modeling of molecular atomization energies with machine learning”  In _Physical Review Letters_ 108.5, 2012, pp. 1–5  DOI: 10.1103/PhysRevLett.108.058301
  * [4] Grégoire Montavon et al.  “Learning Invariant Representations of Molecules for Atomization Energy Prediction”  In _Advances in Neural Information Processing Systems_ 25, 2012, pp. 449–457  URL: https://proceedings.neurips.cc/paper/2012/file/115f89503138416a242f40fb7d7f338e-Paper.pdf
  * [5] Katja Hansen et al.  “Assessment and validation of machine learning methods for predicting molecular atomization energies”  In _Journal of Chemical Theory and Computation_ 9.8, 2013, pp. 3404–3419  DOI: 10.1021/ct400195d
  * [6] Lorenz C. Blum and Jean Louis Reymond  “970 Million druglike small molecules for virtual screening in the chemical universe database GDB-13”  In _Journal of the American Chemical Society_ 131.25, 2009, pp. 8732–8733  DOI: 10.1021/ja902302h
  * [7] Lars Ruddigkeit, Ruud Deursen, Lorenz C. Blum and Jean-Louis Reymond  “Enumeration of 166 Billion Organic Small Molecules in the Chemical Universe Database GDB-17”  In _Journal of Chemical Information and Modeling_ 52.11, 2012, pp. 2864–2875  DOI: 10.1021/ci300415d
  * [8] Felix A. Faber et al.  “Prediction Errors of Molecular Machine Learning Models Lower than Hybrid DFT Error”  In _Journal of Chemical Theory and Computation_ 13.11, 2017, pp. 5255–5264  DOI: 10.1021/acs.jctc.7b00577
  * [9] Nicholas J. Browning, Raghunathan Ramakrishnan, O. Anatole von Lilienfeld and Ursula Roethlisberger  “Genetic Optimization of Training Sets for Improved Machine Learning Models of Molecular Properties”  In _Journal of Physical Chemistry Letters_ 8.7, 2017, pp. 1351–1359  DOI: 10.1021/acs.jpclett.7b00038
  * [10] G. Cybenko  “Approximation by superpositions of a sigmoidal function”  In _Mathematics of Control, Signals and Systems_ 2, 1989, pp. 303–314  DOI: 10.1007/BF02551274
  * [11] Kurt Hornik, Maxwell Stinchcombe and Halbert White  “Multilayer feedforward networks are universal approximators”  In _Neural Networks_ 2.5, 1989, pp. 359–366  DOI: 10.1016/0893-6080(89)90020-8
  * [12] Mohammed N Aldosari, Kiran K Yalamanchi, Xin Gao and S Mani Sarathy  “Predicting entropy and heat capacity of hydrocarbons using machine learning”  In _Energy and AI_ 4 Elsevier Ltd, 2021  DOI: 10.1016/j.egyai.2021.100054
  * [13] Yuanbin Liu, Weixiang Hong and Bingyang Cao  “Machine learning for predicting thermodynamic properties of pure fluids and their mixtures”  In _Energy_ 188, 2019, pp. 116091  DOI: 10.1016/j.energy.2019.116091
  * [14] Matteo Krüger et al.  “Convolutional neural network prediction of molecular properties for aerosol chemistry and health effects”  In _Natural Sciences_ , 2022, pp. e20220016  DOI: https://doi.org/10.1002/ntls.20220016
  * [15] Tobias Fink, Heinz Bruggesser and Jean-Louis Reymond  “Virtual Exploration of the Small-Molecule Chemical Universe below 160 Daltons”  In _Angewandte Chemie International Edition_ 44.10, 2005, pp. 1504–1508  DOI: 10.1002/anie.200462457
  * [16] Tobias Fink and Jean-Louis Reymond  “Virtual Exploration of the Chemical Universe up to 11 Atoms of C, N, O, F: Assembly of 26.4 Million Structures (110.9 Million Stereoisomers) and Analysis for New Ring Systems, Stereochemistry, Physicochemical Properties, Compound Classes, and Drug Discovery”  In _Journal of Chemical Information and Modeling_ 47.2, 2007, pp. 342–353  DOI: 10.1021/ci600423u
  * [17] Daniel Probst and Jean-Louis Reymond  “GDB Databases”, 2022  URL: https://gdb.unibe.ch/downloads/
  * [18] Juntang Zhuang et al.  “AdaBelief Optimizer: Adapting Stepsizes by the Belief in Observed Gradients”  arXiv, 2020  DOI: 10.48550/arXiv.2010.07468
  * [19] Diederik P. Kingma and Jimmy Ba  “Adam: A Method for Stochastic Optimization”  arXiv, 2014  DOI: 10.48550/arXiv.1412.6980
  * [20] Marta Glavatskikh et al.  “Dataset’s chemical diversity limits the generalizability of machine learning predictions”  In _Journal of Cheminformatics_ 11.1, 2019  DOI: 10.1186/s13321-019-0391-2
  * [21] Thomas Cauchy, Benoit Da Mota and Marta Glavastkikh  “PC9_data.zip”, 2019  DOI: 10.6084/m9.figshare.9033977.v1
  * [22] Jeff Bezanson, Alan Edelman, Stefan Karpinski and Viral B. Shah  “Julia: A Fresh Approach to Numerical Computing”  In _SIAM Review_ 59.1, 2017, pp. 65–98  DOI: 10.1137/141000671
  * [23] Mike Innes  “Flux: Elegant machine learning with Julia”  In _Journal of Open Source Software_ 3.25, 2018, pp. 602  DOI: 10.21105/joss.00602
  * [24] Michael Innes et al.  “Fashionable Modelling with Flux”  arXiv, 2018  DOI: 10.48550/arXiv.1811.01457
  * [25] Michael Innes  “Don’t Unroll Adjoint: Differentiating SSA-Form Programs”  arXiv, 2018  DOI: 10.48550/arxiv.1810.07951
  * [26] Mathieu Besançon et al.  “Distributions.jl: Definition and Modeling of Probability Distributions in the JuliaStats Ecosystem”  In _Journal of Statistical Software_ 98.16 Foundation for Open Access Statistic, 2021  DOI: 10.18637/jss.v098.i16
  * [27] Dahua Lin et al.  “JuliaStats/Distributions.jl”  Zenodo, 2022  DOI: 10.5281/zenodo.2647458
  * [28] Tom Breloff  “Plots.jl”  Zenodo, 2022  DOI: 10.5281/zenodo.6989218
  * [29] Tim Holy et al.  “ProgressMeter.jl”, 2021  URL: https://github.com/timholy/ProgressMeter.jl
  * [30] Jarrett Revels et al.  “BenchmarkTools.jl”, 2021  URL: https://github.com/JuliaCI/BenchmarkTools.jl
  * [31] Dahua Lin et al.  “StatsBase.jl”, 2022  URL: https://github.com/JuliaStats/StatsBase.jl
  * [32] Tim Holy et al.  “HDF5.jl”, 2022  URL: https://github.com/JuliaIO/HDF5.jl
  * [33] John M. Kuhn et al.  “Formatting.jl”, 2020  URL: https://github.com/JuliaIO/Formatting.jl
  * [34] Frank R. Burden  “Using Artificial Neural Networks to Predict Biological Activity from Simple Molecular Structural Considerations”  In _Quantitative Structure-Activity Relationships_ 15.1, 1996, pp. 7–11  DOI: https://doi.org/10.1002/qsar.19960150103
  * [35] Frank R. Burden and David A. Winkler  “Robust QSAR Models Using Bayesian Regularized Neural Networks”  In _Journal of Medicinal Chemistry_ 42.16, 1999, pp. 3183–3187  DOI: 10.1021/jm980697n
  * [36] Matthias Rupp  “Machine learning for quantum mechanics in a nutshell”  In _International Journal of Quantum Chemistry_ 115.16, 2015, pp. 1058–1073  DOI: 10.1002/qua.24954
  * [37] L. Riedel  “Ein neues Verfahren zur Abschätzung unbekannter kritischer Drucke von organischen Verbindungen”  In _Zeitschrift für Elektrochemie und angewandte physikalische Chemie_ 53.4, 1949, pp. 222–228  DOI: 10.1002/bbpc.19490530411
  * [38] Aksel Lydersen  “Estimation of Critical Properties of Organic Compounds”  In _Engineering Experimental Station Report 3_ 3, 1955
  * [39] D. Ambrose  “Correlation and Estimation of Vapor-Liquid Critical Properties: I. Critical Temperatures of Organic Compounds”, 1978, corrected 1980
  * [40] D. Ambrose  “Correlation and Estimation of Vapor-Liquid Critical Properties: II. Critical Pressures and Volumes of Organic Compounds”, 1979
  * [41] K. G. Joback and R. C. Reid  “Estimation of Pure-Component Properties From Group-Contributions”  In _Chemical Engineering Communications_ 57.1-6 Taylor & Francis, 1987, pp. 233–243  DOI: 10.1080/00986448708960487
  * [42] Leonidas Constantinou and Rafiqul Gani  “New group contribution method for estimating properties of pure compounds”  In _AIChE Journal_ 40.10, 1994, pp. 1697–1710  DOI: 10.1002/aic.690401011
  * [43] Yash Nannoolal, Jürgen Rarey and Deresh Ramjugernath  “Estimation of pure component properties: Part 2. Estimation of critical property data by group contribution”  In _Fluid Phase Equilibria_ 252.1, 2007, pp. 1–27  DOI: 10.1016/j.fluid.2006.11.014
  * [44] Yash Nannoolal, Jürgen Rarey, Deresh Ramjugernath and Wilfried Cordes  “Estimation of pure component properties: Part 1. Estimation of the normal boiling point of non-electrolyte organic compounds via group contributions and group interactions”  In _Fluid Phase Equilibria_ 226, 2004, pp. 45–63  DOI: 10.1016/j.fluid.2004.09.001
  * [45] Sidney W. Benson and Jerry H. Buss  “Additivity Rules for the Estimation of Molecular Properties. Thermodynamic Properties”  In _The Journal of Chemical Physics_ 29.3, 1958, pp. 546–572  DOI: 10.1063/1.1744539
  * [46] Aage Fredenslund, Russell L. Jones and John M. Prausnitz  “Group-contribution estimation of activity coefficients in nonideal liquid mixtures”  In _AIChE Journal_ 21.6, 1975, pp. 1086–1099  DOI: 10.1002/aic.690210607
  * [47] Guido Van Rossum and Fred L Drake Jr  “Python tutorial”  Centrum voor Wiskunde en Informatica Amsterdam, The Netherlands, 1995
  * [48] Anders S. Christensen et al.  “qmlcode/qml: Release v0.3.1”  Zenodo, 2017  DOI: 10.5281/zenodo.817332
  * [49] Travis E. Oliphant  “Python for Scientific Computing”  In _Computing in Science & Engineering_ 9.3, 2007, pp. 10–20  DOI: 10.1109/MCSE.2007.58
  * [50] Pavlo O. Dral, O. Anatole von Lilienfeld and Walter Thiel  “Machine learning of parameters for accurate semiempirical quantum chemical calculations”  In _Journal of Chemical Theory and Computation_ 11.5, 2015, pp. 2120–2125  DOI: 10.1021/acs.jctc.5b00141
  * [51] Sidney W. Benson  “Bond Energies”  In _Journal of Chemical Education_ 42.9, 1965, pp. 502–518  DOI: 10.1021/ed042p502
  * [52] James J. P. Stewart  “Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements”  In _Journal of Molecular Modeling_ 13.12, 2007, pp. 1173–1213  DOI: 10.1007/s00894-007-0233-4
  * [53] T. M. Mitchell  “Machine Learning”, McGraw-Hill International Editions  McGraw-Hill, 1997
  * [54] Xavier Glorot and Yoshua Bengio  “Understanding the difficulty of training deep feedforward neural networks”  In _Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics_ 9, Proceedings of Machine Learning Research  Chia Laguna Resort, Sardinia, Italy: PMLR, 2010, pp. 249–256  URL: https://proceedings.mlr.press/v9/glorot10a.html
  * [55] Sebastian Ruder  “An overview of gradient descent optimization algorithms”, 2016, pp. 1–14  arXiv: https://arxiv.org/abs/1609.04747
  * [56] John Duchi, Elad Hazan and Yoram Singer  “Adaptive Subgradient Methods for Online Learning and Stochastic Optimization”  In _Journal of Machine Learning Research_ 12, 2011, pp. 2121–2159
  * [57] Alex Graves  “Generating Sequences With Recurrent Neural Networks”  arXiv, 2013  DOI: 10.48550/arXiv.1308.0850
  * [58] Ashia C. Wilson et al.  “The Marginal Value of Adaptive Gradient Methods in Machine Learning”  arXiv, 2017  DOI: 10.48550/arXiv.1705.08292
  * [59] Thomas Elsken, Jan Hendrik Metzen and Frank Hutter  “Neural Architecture Search: A Survey”  In _Journal of Machine Learning Research_ 20.55, 2019, pp. 1–21  URL: http://jmlr.org/papers/v20/18-598.html
  * [60] Martin Wistuba, Ambrish Rawat and Tejaswini Pedapati  “A Survey on Neural Architecture Search”  arXiv, 2019  DOI: 10.48550/arXiv.1905.01392
  * [61] George S. Hammond  “A Correlation of Reaction Rates”  In _Journal of the American Chemical Society_ 77.2, 1955, pp. 334–338  DOI: 10.1021/ja01607a027
  * [62] A. Einstein  “Die Plancksche Theorie der Strahlung und die Theorie der spezifischen Wärme”  In _Annalen der Physik_ 327.1, 1907, pp. 180–190  DOI: 10.1002/andp.19063270110
  * [63] P. Debye  “Zur Theorie der spezifischen Wärmen”  In _Annalen der Physik_ 344.14, 1912, pp. 789–839  DOI: 10.1002/andp.19123441404
  * [64] Katja Hansen et al.  “Machine Learning Predictions of Molecular Properties: Accurate Many-Body Potentials and Nonlocality in Chemical Space”  In _The Journal of Physical Chemistry Letters_ 6.12, 2015, pp. 2326–2331  DOI: 10.1021/acs.jpclett.5b00831
  * [65] Bing Huang and O. Anatole von Lilienfeld  “Communication: Understanding molecular representations in machine learning: The role of uniqueness and target similarity”  In _The Journal of Chemical Physics_ 145.16, 2016, pp. 161102  DOI: 10.1063/1.4964627
  * [66] Jörg Behler and Michele Parrinello  “Generalized neural-network representation of high-dimensional potential-energy surfaces”  In _Physical Review Letters_ 98.14, 2007, pp. 1–4  DOI: 10.1103/PhysRevLett.98.146401
  * [67] Jörg Behler  “Atom-centered symmetry functions for constructing high-dimensional neural network potentials”  In _Journal of Chemical Physics_ 134.7, 2011  DOI: 10.1063/1.3553717
  * [68] J. S. Smith, O. Isayev and A. E. Roitberg  “ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost”  In _Chemical Science_ 8.4 Royal Society of Chemistry, 2017, pp. 3192–3203  DOI: 10.1039/C6SC05720A