a an the of on at with for and in into from by without DFG Cluster of
Excellence “Unifying Systems in Catalysis” (UniSysCat), Technische Universität
Berlin, 10623 Berlin, Germany Department of Chemistry, Brown University,
Providence RI, USA

# Impact of the characteristics of quantum chemical databases on machine
learning predictions of tautomerization energies.

Luis Itza Vazquez-Salazar  Department of Chemistry, University of Basel,
Klingelbergstrasse 80 , CH-4056 Basel, Switzerland  Eric Boittier  Department
of Chemistry, University of Basel, Klingelbergstrasse 80 , CH-4056 Basel,
Switzerland  Oliver T. Unke  Machine Learning Group, Technische Universität
Berlin, 10587 Berlin, Germany  Markus Meuwly  Department of Chemistry,
University of Basel, Klingelbergstrasse 80 , CH-4056 Basel, Switzerland.
<EMAIL_ADDRESS>
###### Abstract

An essential aspect for adequate predictions of chemical properties by machine
learning models is the database used for training them. However, studies that
analyze how the content and structure of the databases used for training
impact the prediction quality are scarce. In this work, we analyze and
quantify the relationships learned by a machine learning model (Neural
Network) trained on five different reference databases (QM9, PC9, ANI-1E,
ANI-1 and ANI-1x) to predict tautomerization energies from molecules in
Tautobase. For this, characteristics such as the number of heavy atoms in a
molecule, number of atoms of a given element, bond composition, or initial
geometry on the quality of the predictions are considered. The results
indicate that training on a chemically diverse database is crucial for
obtaining good results but also that conformational sampling can partly
compensate for limited coverage of chemical diversity. The overall best
performing reference database (ANI-1x) performs on average by 1 kcal/mol
better than PC9 which, however, contains about two orders of magnitude fewer
reference structures. On the other hand, PC9 is chemically more diverse by a
factor of $\sim 5$ as quantified by the number of amons it contains compared
with the ANI family of databases. We explicitly demonstrate that when certain
types of bonds need to be covered in the target database (Tautobase) but are
undersampled in the reference databases the resulting predictions are poor.
Examples include C(sp2)-C(sp2) double bonds close to hetero atoms and azoles
containing N-N and N-O bonds. A quantitative measure for these deficiencies is
the Kullback-Leibler divergence between reference and target distributions.
Analysis of the results with a TreeMAP algorithm provides deeper understanding
of specific deficiencies in the reference data sets. Capitalizing on this
information can be used to either improve existing databases or to generate
new databases of sufficient diversity for a range of ML applications in
chemistry.

##  1 Introduction

In the last decade, the application of machine learning (ML) techniques in
chemistry has significantly increased1, 2, 3, 4. This has occasionally been
related to a paradigm shift, revolutionizing the available techniques to
understand and simulate chemistry5, 6. The excitement is seemingly justified,
given the outcomes of ML techniques’ central promise that, by using a
sufficiently large number of examples and a rule-discovery algorithm, it is
possible to obtain a scientific understanding of the underlying relationships
covered by the data3, 1. Furthermore, ML techniques are fast compared with
quantum chemical methods, while also reaching comparable accuracy7, 8, 9, 10,
11, 12, 13, 14, 15, 16, 17.

On the other hand, application of quantum ML methods to concrete problems
requires large amounts of data which first need to be generated from
electronic structure calculations18, 19, 20. Consequently, data generation is
computationally demanding. An essential challenge for the extension of ML
methods’ applicability in chemistry is understanding how suitable databases
can be constructed to maximize accuracy and transferability of the models. An
important ingredient for this step is the degree and confidence with which a
human can understand the relationship between cause (starting database and
model) and result or observation (applying the model to a new task)21, 22.
This process has also been called “interpretability” and it can be used to
understand the relationships learned by the model or contained in the data
used for training it23, 24. Part of the present work is concerned with the aim
to relate the composition of the initial chemical databases based on which ML
models are conceived with their performance on the prediction of a property of
interest (tautomerization energy) on a set of unseen examples.

To test the effect of different databases on the reliability of the ML model,
the problem of predicting tautomerization energies is considered. Tautomerism
is a form of reversible isomerization involving the rearrangement of a charged
leaving group within a molecule25 (e.g. Fig. 1a). One isomer transforms into
the other by a heterolytic splitting followed by a recombination of the
fragments formed26. This process involves the migration of one or more double
bonds and atoms or groups. The isomers (i.e. tautomers) generated in this
reaction are chemically independent species with defined properties27. It is
known that this type of reaction is of importance for biological molecules
such as amino acids26, DNA28, 29, RNA30, and atmospheric processes.31
Additionally, it is estimated that tautomerism can occur in up to two thirds
of small molecules32, and a majority of commercial drugs33, 34.

Despite its widespread occurrence and importance, quantitative studies of
tautomerism are still challenging because small changes in molecular structure
or solvent environment can dramatically change the tautomeric equilibrium27,
35. Moreover, small free energy differences between two tautomers in solution
make the use of high level theoretical methods and an adequate basis set
mandatory which limits its use for calculations of tautomerization energies
and ratios. 35, 36 As an example, tautomerization in malonaldehyde (MA) is
considered. MA has served as a prime example to develop and test computational
methods for a realistic description of hydrogen transfer in small molecules.37
Experimentally, the ground state tunneling splitting is 21.58314 cm-1 which
has been determined by different experiments with very high accuracy. 38, 39
Furthermore, proton transfer rates in a di-imine derivative have been
determined with nuclear magnetic resonance (NMR) spectroscopy. 40 Such
experiments provide direct information on the barrier height separating the
two tautomeric states “A” and “B”. Using a state-of-the art full-dimensional
potential energy surface at the near basis-set-limit frozen-core CCSD(T) level
of theory,41 the tunneling splitting from quantum simulations was determined
as 23.4 cm-1.42 Alternatively, using a reduced dimensionality Hamiltonian, the
barrier height for proton transfer in a parametrized molecular mechanics with
proton transfer (MMPT) potential was found to be 4.34 kcal/mol which yields a
tunneling splitting of 21.2 cm-1, consistent with experiment.43, 44 This
barrier height is close to the value from CCSD(T) calculations which yield 4.1
kcal/mol.41 These examples illustrate that calculations at the highest levels
of theory are required for quantitative studies of the energetics underlying
tautomerization.

In the last decade, development of ML models has allowed the design of robust
models that can routinely reach prediction errors better than chemical
accuracy at low computational cost12, 45. However, there have been few
discussions on how databases can be improved/designed to obtain better
predictions from the ML model. Ideally, the combination of a robust ML model
and an adequate database will result in quantitative results for the
prediction of a property of interest. The availability of public databases of
tautomers46, 47 makes the prediction of tautomerization energies using these
ML models an ideal test case to study how different training databases
influence the accuracy of ML methods.

The present work is structured as follows. First, the methods, databases, and
the analysis performed are introduced. Next, the results for the
tautomerization energy predictions using models trained on the different
tested databases are presented. Additionally, prediction errors for
tautomerization energies are analyzed. The effect of different characteristics
of the training data on predicting the tautomerization energy and the
individual molecules’ energy are evaluated. Finally, the results are discussed
and conclusions regarding the findings and interpretability of broadly
conceived and learned ML models applied to a specific chemical question are
drawn.

Figure 1: (A) Tautomerism is a form of reversible isomerization involving the
rearrangement of a charged leaving group within a molecule. The keto-enol
tautomerism of acetone, an equilibrium which heavily favours the keto side, is
shown as an example. (B) Chemical space can be decomposed systematically
through the use of atoms-in-molecules-based fragments (amons).48 The amons
present in acetone (SMILES: CC(=O)C), as well as their corresponding SMILES
are given as an example.

##  2 Methods

###  2.1 Machine Learning

PhysNet was used for the representation and evaluation of the data sets, using
the hyperparameters from the original publication45. For training, only the
nuclear charges ($Z$), the energies of the molecules ($E$) and their
coordinates ($\mathbf{R}$) were considered. The energies used for training
were those reported by the different databases minus the atomization energy at
the given level of theory. In all cases, a training, testing and validation
split of 8 : 1 : 1 was used. The loss function was

$\mathcal{L}=w_{\rm{E}}|E-E^{\rm{ref}}|+\lambda_{\rm{nh}}\mathcal{L}_{\rm{nh}}$
(1)

where $E^{\rm{ref}}$ is the reference energy, $w_{\rm{E}}=1$ is the weighting
hyperparameter for the energy, $\lambda_{\rm{nh}}=10^{-2}$ is a regularization
hyperparameter and the term $\mathcal{L}_{\rm{nh}}$ is a non hierachicality
regularization penalty. The loss function (Eq. 1) was minimized using AMSgrad
with a learning rate of $10^{-3}$. Overfitting was prevented by the use of
early stopping, the convergence criteria considered is the saturation of the
validation-loss function45.

###  2.2 Database Selection

For training the NNs, four widely used databases for benchmarking predictive
models of DFT-based energies are employed, namely QM949, PC950, ANI-151 and
ANI-1x.52 An additional database, ANI-1E (where "E" stands for equilibrium)
containing only the equilibrium structures of the ANI-1 database, was
generated at the $\omega$B97x/6-31G(d) level of theory. These databases can be
divided into categories based on the type of geometries they contain. The
datasets consisting solely of equilibrium structures are QM9, PC9 and ANI-1E,
and sample only chemical space. Contrary to that, the ANI-1 and ANI-1x
databases contain equilibrium and non-equilibrium structures which sample
chemical and conformational space.

Training sets: The QM9 data set49 was generated as a subset of the GDB-17
chemical universe53, consisting of 133885 molecules, containing less than or
equal to 9 heavy atoms (either C, N, O or F). Reference energies were computed
at the B3LYP/6-31G(2df,p) level of theory. For the present work, QM9 was
filtered to include only molecules which passed a geometry consistence check49
and considering only those containing carbon, nitrogen or oxygen atoms. The
final size of the QM9 training dataset used here consisted of 110426
molecules.

The PC9 dataset was created as an alternative to QM9 to improve coverage of
chemical space.50 It is a subset of the PubChemQC54 and is limited to
molecules with 9 heavy atoms or less ($n_{\rm atoms}\leq 9$). This database
consists of 99234 molecules, calculated at the B3LYP/6-31G(d) level of theory,
and excludes enantiomers, tautomers, isotopes as well as other specific
artifacts in PubChemQC50. PC9 also contains 5325 molecules with an electronic
state different from a singlet which were removed for the present work. As in
the case of QM9, molecules which contain fluorine were removed. The final size
of this dataset was 85875 molecules.

ANI-151 consists of 24 million geometries generated using normal mode sampling
from 57462 unique molecules. ANI-1 is a subset of the GDB-11 chemical
universe55, 56. A related dataset, ANI-1x 52, was created using an active
learning57 procedure which reduced the original ANI-1 database to 5 million
structures. Starting from the ANI-1 database51, 58 the ANI-1E dataset was
generated and consists only of the corresponding equilibrium structures. The
new ANI-1E database contains 57462 molecules limited to eight heavy atoms
(either C, O or N). The generation of this database is further described below
in the subsection of electronic structure calculations.

Tautomerization energy evaluation set: The performance of the NN models
described above was evaluated on a subset of molecules from TautoBase46, a
public database of 1680 tautomer pairs. The Tautobase was filtered to
molecules only containing hydrogen, carbon, nitrogen, or oxygen atoms. The
size of the final test set was 1257 tautomer pairs (2514 molecules). The
geometry generation and structural optimization for these molecules is
described below.

###  2.3 Initial geometry

To investigate the effect of the geometry of the molecules passed to the NN
model on its performance, a second set of geometries for the tautobase was
also evaluated. These geometries were generated from the SMILES representation
using OpenBabel59 and were optimised with the MMFF94 force field60.

Additionally, a subset of the test set composed of 34 tautomeric pairs which
were part of the SAMPL2 challenge61 were considered. Those 34 pairs (68
molecules) were optimized by six popular general atomistic force fields:
CHARMM2762, GAFF63, OPLS64, UFF65, Gromos66, and Ghemical67. Details on the
generation of the geometries are reported in the SI.

###  2.4 Electronic Structure Calculations

Generation of ANI-1E: Starting from the ANI-1 database58, a new data set,
ANI-1E, was generated. From the SMILES strings provided by Smith et al. 58
initial geometries using OpenBabel59 were generated. Subsequently, geometries
were optimised using PM768 implemented in MOPAC201669, before a final geometry
optimization and frequency calculation at the $\omega$B97x/6-31G(d) level of
theory was performed using Gaussian0970. The final results were checked to
assure they did not contain imaginary frequencies and therefore correspond to
minima on the potential energy surface.

Tautomerization evaluation set: The molecules used for the evaluation of the
NN models were generated from the SMILES provided in Ref. 46 using the
OpenBabel software59. These structures were then optimized at the same level
of theory for each of the databases used for training (QM9:
B3LYP/6-31G(2df,p), PC9: B3LYP/6-31G(d), ANI: $\omega$B97x/6-31G(d)) using
Gaussian0970. Here, the tautomerization energy $\Delta E_{\rm Tauto}$ is
defined as the energy difference between tautomers A and B in their optimized
structures, see Figure 1. These optimized geometries were given as input to
the respective, previously trained NN models using the five reference
databases QM9, PC9, ANI-1E, ANI-1 and ANI-1x.

Table 1: Overview of the training datasets used in this work. QM9, PC9 and ANI-1E contain equilibrium structures and can be considered to only sample chemical space, whereas ANI-1 and ANI-1x also contain non-equilibrium geometries which sample conformational space. The number of molecules refers to the total number of data in each dataset.aStructures generated through normal mode sampling. bTraining set selected using active learning. Database | Number of Molecules | Level of Theory | Parent Universe |
---|---|---|---|---
QM9 | 128908 | B3LYP/6-31G(2df,p) | GDB-17 |
PC9 | 85870 | B3LYP/6-31G(d) | PubChemQC |
ANI-1E | 57462 |  $\omega$B97x/6-31G(d) | GDB-11 |
ANI-1a | 24 million |  $\omega$B97x/6-31G(d) | GDB-11 |
ANI-1xb | 5 million |  $\omega$B97x/6-31G(d) | ANI-1 |

###  2.5 Comparison of structural properties of different databases

As a way to compare the composition of the different datasets evaluated in
terms of structural properties (e.g. bond lengths), a Gaussian kernel density
estimation71 of their distributions was generated, see Figures
LABEL:sifig:dist_bonds_C_9 to LABEL:sifig:dist_bonds_N_9p. The similarities
between the distributions used to train the NN models and those from the test
set of tautomers was quantified by computing the relative entropy (or
Kullback-Leibler (KL) divergence)72

$D(p\parallel
q)=\int_{-\infty}^{\infty}p(x)log\left(\frac{p(x)}{q(x)}\right)dx$ (2)

This metric quantifies the overlap between a reference distribution $p(x)$ and
a target distribution $q(x)$. Because the KL divergence is not symmetric
($D(p\parallel q)\neq D(q\parallel p)$), it is important to specify which
distribution is used as the reference. In the present work, the Tautobase
(target) distribution and the QM9, PC9, ANI-1E are the reference
distributions. The KL divergence allows to quantify how much information of
the reference databases (i.e. QM9, PC9, ANI-1E) is ’missing’ to best cover the
information contained in Tautobase.

###  2.6 Chemical Space “Coverage” from Fragment Analysis

The ‘coverage’ of chemical space contained in the Tautobase, QM9, PC9 and ANI
family of databases was analysed with respect to the atom-in-molecule-based
fragments (amons)48. The amons are generated from the SMILES representation of
the molecule, see Figure 1b. This representation is used to construct a
molecular graph from which sub-graphs to a maximum number of atoms (excluding
hydrogen) are generated. All sub-graphs are checked to be valid and unique.
Here, amons up to and including a maximum of five heavy atoms were generated
by an in-house script.

##  3 Performance of the NN on the Tautobase

###  3.1 Overall Performance

The mean absolute errors for the tautomerization energies $\Delta E_{\rm
Tauto}$ range from 1.68 kcal/mol (ANI-1x) to 4.59 kcal/mol (ANI-1). The
results are summarized in Table 2 for all molecules in the test set and
graphically reported in Figures 2 (datasets with equilibrium structures) and 3
(datasets with both equilibrium and non-equilibrium structures). The
prediction errors for the energy of single isomers, $E_{\rm SI}$, with respect
to DFT energies are also reported (Table 2). Note that because the
tautomerization energy is defined as the difference between the isomer
energies, predictions of tautomerization energies are often more accurate due
to cancellation of systematic errors. On the other hand, the energies for
single isomers are considerably larger and span a much wider range because
they scale with the number of atoms that make up a molecule. Therefore, the
NN-based energies for larger molecules are expected to be associated with
considerably larger errors.

Figure 2: Correlation between the calculated (DFT) and predicted (NN)
tautomerization energies ($\Delta E_{\rm{Tauto}}=E_{\textrm{A}}-E_{\rm{B}}$)
for molecules with $n_{\rm{atoms}}\leq 9$ (left) and $n_{\rm{atoms}}>9$
(right) for the models trained on the datasets which only cover chemical
space. Pearson correlation coefficients ($r^{2}$) and the Mean Absolute Error
(MAE) values are reported in brackets as [$r^{2}$, MAE].

To assess whether the accuracy for predicting tautomerization energies
correlates with the performance of the trained NNs on the chemical databases,
the QM9, PC9, and ANI-1E models are considered. The MAEs on held-out test sets
for each of the respective training runs are 0.10 kcal/mol (QM9), 0.69
kcal/mol (PC9), and 0.27 kcal/mol (ANI-1E) which is comparable to values
ranging from 0.19 kcal/mol to 0.30 kcal/mol for the QM9 data set (depending on
the sizes of the training and validation sets used).45 However, when applying
these trained NNs to evaluate Tautobase, the MAEs are 3.67 kcal/mol (QM9),
2.60 kcal/mol (PC9), and 3.66 kcal/mol (ANI-1E), respectively. Hence, there
appears to be no correlation between the quality of the trained NNs (measured
on test data sampled from the same database used for training) and their
performance on Tautobase.

Table 2: Mean Absolute (MAE) and Root-Mean-Squared Error (RMSE) for the prediction of tautomerization energy $\Delta E_{\rm Tauto}$, and the single isomer energies, $E_{\rm SI}$, for the entire Tautobase (1257 tautomeric pairs) for each of the datasets. | $\Delta E_{\textrm{Tauto}}$ | $E_{\textrm{SI}}$
---|---|---
Database | MAE | RMSE | MAE | RMSE
QM9 | 3.67 | 7.12 | 5.00 | 8.40
PC9 | 2.60 | 5.41 | 6.90 | 13.20
ANI-1E | 3.66 | 7.09 | 15.20 | 17.50
ANI-1 | 4.59 | 7.56 | 13.40 | 17.00
ANI-1x | 1.68 | 2.85 | 1.80 | 3.60

Next, the performance of the trained models for predicting $\Delta E_{\rm
Tauto}$ and $E_{\rm SI}$ depending on the number of heavy atoms is assessed.
For this, results for the subset of molecules with $n_{\rm atoms}\leq 9$,
referred to as “Set1” in the following, is considered separately from those
with $n_{\rm atoms}>9$, which is “Set2”. This distinction is motivated by the
fact that the PC9 and QM9 databases contain structures with only up to 9 heavy
atoms, i.e. models need to extrapolate for larger structures. For Set1, the
PC9 (Figure 2C) and ANI-1x (Figure 3C) data sets perform best. Both achieve
chemical accuracy (MAE < 1 kcal/mol) with respect to the DFT values for
$\Delta E_{\rm Tauto}$.

Figure 3: Correlation between the calculated (DFT) and predicted (NN)
tautomerization energies ($\Delta E_{\rm{Tauto}}=E_{\textrm{A}}-E_{\rm{B}}$)
for molecules with $n_{\rm{atoms}}\leq 9$ (left) and $n_{\rm{atoms}}>9$
(right) for the models trained with the datasets which cover chemical and
conformational space. Pearson correlation coefficients ($r^{2}$) and Mean
Absolute Error (MAE) values are reported in brackets [$r^{2}$, MAE].

The extrapolation to Set2 increases the prediction errors for most of the
databases studied. Again, the ANI-1x database performed best for $\Delta
E_{\rm{Tauto}}$ with a MAE of 2.20 kcal/mol, followed by PC9, QM9, ANI-1E and
ANI-1 with a MAE of 6.29 kcal/mol. The number of atoms in the database also
influences the correlation coefficient $r^{2}$. A better correlation is
observed when the size of evaluated molecules is in the range covered by the
training database, i.e. for Set1. Here, the correlation coefficients range
from 0.69 (QM9) to 0.99 (ANI-1x). For Set2, the $r^{2}$ values are
significantly lower. A particularly noteworthy case is the QM9 database, which
shows almost the same MAE and $r^{2}$ for both subsets of the tautobase
(Figures 2A and B). The performance for Set1 and Set2 also differs in the
number and magnitude of outliers, see Figures 2 and 3.

The RMSE in Table 2 show that the spread for $E_{\rm SI}$ could be a reason
for large outliers (see Figure LABEL:sifig:error_dist). It is likely that
there is some error cancellation when predicting tautomerization energies
(i.e. energy differences). For example, if the trained NN predicts too large
energies for both isomers, this systematic error cancels when their energy
difference is computed.

###  3.2 Error analysis

Next, prediction errors for $\Delta E_{\rm Tauto}$ and $E_{\rm SI}$ are
analyzed and discussed for all trained models. In Figures 4A and B, the kernel
density estimate of the error distribution is reported. The violin plots in
Figures 4C and D show the spread of errors, which helps to identify large
outliers.

Figure 4: Error analysis on the prediction of the tautomerization energies.
Panels A and B: Kernel density estimate of the error on the prediction of the
tautomerization energies for the different databases evaluated in this work.
Panels C and D: Normalized error distribution up to the 95% quantile of the
different datasets for the errors in the tautomerization energy. The blackbox
inside spans between the 25% and 75% quantiles with a white dot indicating the
mean of the distribution. The whiskers indicate the 5% and 95 % quantiles. The
panels on the left and right are for Set1 and Set2, respectively.

In all cases studied, the distribution of errors for $\Delta\Delta E_{\rm
Tauto}=\Delta E_{\rm Tauto}^{\rm DFT}-\Delta E_{\rm Tauto}^{\rm NN}$ is
centered around zero. The width of the distribution depends on the reference
dataset and on the number of atoms in the molecule. For Set1 (Figure 4A), the
error distributions for PC9 and ANI-1x suggest high probabilities
($p(\Delta\Delta E_{\rm Tauto})=p(\Delta E_{\rm Tauto}^{\rm DFT}-\Delta E_{\rm
Tauto}^{\rm NN})$), around 60 %, to obtain a small error. Conversely, QM9
performs worst with a maximum height of only $p(\Delta\Delta E_{\rm
Tauto})\sim 15$ % and a faint maximum below $\Delta\Delta E_{\rm Tauto}=-10$
kcal/mol to predict such energy differences with a larger probability than a
positive value. The error distributions for ANI-1E and ANI-1 are similar in
shape which indicates that their performance is comparable although the number
of structures in ANI-1E is one order of magnitude smaller than that in ANI-1.
Hence, adding additional structures (ANI-1 vs. ANI-1E) does not necessarily
improve performance.

On the other hand, for Set2 (Figure 4B), the performance of QM9 and PC9 is
comparable given the similar shape of their distribution of $p(\Delta\Delta
E_{\rm Tauto})$, ANI-1x gives the best predictions with a maximum height of
around 15 % to obtain an error for $\Delta E_{\rm Tauto}$ close to zero. All
other reference data sets perform inferior with ANI-1 reaching only a 5 %
$p(\Delta\Delta E_{\rm Tauto})$ for a prediction close to zero. In addition,
for most of the data sets the error distribution is asymmetric with an
increased probability to predict a negative value for $\Delta\Delta E_{\rm
Tauto}$ compared to a positive value.

Results for the normalized error distributions $p(\Delta\Delta E_{\rm Tauto})$
are shown in Figures 4C and D. For Set1, PC9 and ANI-1x show the smallest
outliers by magnitude with an error below 2.5 kcal/mol. On the other hand, QM9
has the largest outliers with some errors larger than 15 kcal/mol. The average
error for all reference distributions is around or below 1 kcal/mol for the 75
% quantile. For molecules in Set2, ANI-1 has the largest outliers, followed by
ANI-1E, QM9, PC9, and ANI-1x performing best with a maximum error of around 5
kcal/mol, see Figure 4D.

For completeness, error distributions $p(\Delta E_{\rm SI})=p(E_{\rm SI}^{\rm
DFT}-E_{\rm SI}^{\rm NN})$ for individual molecules and their normalized
variants are also reported in Figures LABEL:sifig:error_distA to D. For Set1,
the distributions for PC9 and ANI-1x are centered around zero with peak
heights at 80 % which decreases to 25 % for ANI-1E. For ANI-1 it is shifted to
negative and for QM9 to positive values. For Set2 (Figure
LABEL:sifig:error_distB), all error distributions are asymmetric and extend to
large negative values of $\Delta E_{\rm SI}$. The best and worst performing
reference distributions are ANI-1x and ANI-1, respectively. The normalized
error distributions (Figures LABEL:sifig:error_distC and D) for both sets are
strongly peaked. For Set1 (Figure LABEL:sifig:error_distC) the maxima for PC9
(2.5 kcal/mol) is the lowest whereas ANI-1 has the largest errors. For Set2
(Figure LABEL:sifig:error_distD), the outliers are even more pronounced with
$|\Delta E_{\rm SI}|>100$ kcal/mol for ANI-1E and ANI-1. In general, the
performance of ANI-1E is better than that of ANI-1 with a smaller MAE,
outliers of smaller magnitude and a more compact distribution. These results
are surprising, given the large difference between the size of both datasets
(ANI-1E ($\approx 57$k) and ANI-1 ($\approx 20$M)) and confirm the earlier
observation that addition of new structures to a database does not necessarily
improve performance.

In summary, for Set1 the database with broader chemical diversity (PC9) and
the database with the widest sampling of chemical and conformational space
(ANI-1x) perform best. Hence, chemical diversity is essential for faithful
prediction of $\Delta E_{\rm Tauto}$ but it can be substituted to some extent
with adequate sampling of conformational space. For larger molecules (Set2),
the best results are obtained by ANI-1x which suggests that sampling of
conformational space improves extrapolation to larger molecules. Datasets
containing only equilibrium structures perform similarly for predicting
$\Delta E_{\rm Tauto}$.

##  4 Effect of different database characteristics on predictions

This section analyzes the predictive power of the NNs trained on the five
different training databases for tautomerization energies by considering
various chemical properties such as the number of heavy atoms, the number of
atoms of a given element, or the type of chemical bonds. Given the non-linear
nature of the NN, the relationships between these characteristics and
whether/how they are related to the performance of the model is a challenging
task. The features studied here were selected because they might be considered
for the selection of a training database for the prediction of a chemical
property (in this case the tautomerization energy) or because they can be
optimized for the generation/enhancement of datasets used to train models for
specific purposes.

###  4.1 Number of atoms

The first characteristic considered was the number of heavy atoms (C, N and O)
contained in the reference data sets and how this affects the prediction
quality on Tautobase. For Set1 the MAE for $\Delta E_{\rm Tauto}$ typically
decreases with increasing molecular size, see Figure 5A, for all five
reference data sets. This can be broadly related to the increase in the number
of molecules with the number of heavy atoms contained in the reference
databases used for training (see Figure LABEL:sifig:NatomsperDatabases). For
all data sets except for QM9, the MAE decreases to levels below 0.5 kcal/mol
as the number of heavy atoms increases.

For larger molecules (Set2) in the Tautobase, the MAEs increase significantly,
see Figure 5B. Broadly speaking, for up to 25 heavy atoms the MAE is still
within 5 kcal/mol but increases considerably for larger molecules. ANI-1x
performs best with ${\rm MAE}<1$ kcal/mol up to $n_{\rm atoms}=25$ but errors
increase above 10 kcal/mol beyond that. This is followed by QM9 and PC9 which,
on average, have MAEs of $\sim 2$ kcal/mol followed by ANI-1E and ANI-1.

Figure 5: Mean Absolute Error (MAE) by number of heavy atoms (C,O,N) on the
molecule for the tautomerization energy. The number of molecules for
increasing number of heavy atoms is shown as a histogram. Panel A: Results for
Set1, i.e. molecules with $n_{\rm{atoms}}\leq 9$. Panel B: for Set2. The inset
in panel B shows the MAE for $\Delta E_{\rm Tauto}$ for $10\leq n_{\rm
atoms}\leq 25$.

The MAE for predicting $E_{\rm SI}$ for Set1 (Figure
LABEL:sifig:MAEvNAtom_molA) shows a slightly different profile than for
$\Delta E_{\rm Tauto}$ in that databases containing only equilibrium
structures show large errors for molecules with 3 and 4 heavy atoms. With
increasing size the error decreases. This is most pronounced for PC9 which
eventually achieves the same quality as ANI-1x. It should be noted that for
$E_{\rm SI}$ the overall shape of the profile of MAE vs. $n_{\rm atoms}$ for
databases which only contain equilibrium structures is similar but the
magnitude of the MAE differs. This is a consequence of the chemical diversity
of the databases as discussed in subsection 4.4. For ANI-1 and Set1, the MAE
is smallest for $n_{\rm atoms}=5$ and then starts to grow again. For Set2
(Figure LABEL:sifig:MAEvNAtom_molB), the MAE displays a steady increase with
the number of heavy atoms.

In summary, ANI-1x performs best across all values of $n_{\rm atoms}$ for the
Tautobase, followed by PC9 across most values for $n_{\rm atoms}$. For Set2,
QM9 is quite reliable whereas ANI-1 and ANI-1E perform worst which reiterates
the earlier finding that adding perturbed structures to a data set does not
necessarily improve the quality on the task at hand (which is the estimation
of $\Delta E_{\rm Tauto}$). Consequently, the results can be worse than those
obtained when training only on equilibrium structures.

A further refinement can be made by analyzing the predictions for $\Delta
E_{\rm Tauto}$ in terms of the number of C-, N-, and O-atoms contained in the
molecules of the reference database (Tautobase), see Figure 6. For Set1 the
prediction error tends to decrease (except for QM9) with increasing number of
carbon atoms as shown in Figure 6A whereas for Set2 it increases to different
extents depending on the reference database considered, see Figure 6B. For
nitrogen and oxygen atoms and ANI-1x all MAEs for Set1 and Set2 are small
($\sim 1$ kcal/mol), except for the largest numbers of N-atoms, see Figure 6F.
For the PC9, ANI-1E, and ANI-1 databases and Set1 all MAEs are below or around
1 kcal/mol whereas for QM9 they can be larger. For Set2, the MAEs are up to 5
kcal/mol for molecules for which at least tens of representatives are
contained in Tautobase, but start to increase significantly below that, see
Figure 6D and F.

Figure 6: Mean Absolute Error (MAE) by number of atoms of a given element for
the tautomerization energy. A histogram of the number of molecules for
different numbers of heavy atoms is shown in the background. Panels A and B
show the results by number of carbon atoms. Panels C and D shows the results
by number of nitrogen atoms. Finally, panels E and F show the results by
number of oxygen atoms. Left panel shows results for molecules with
$n_{\rm{atoms}}\leq 9$ the right for molecules with $n_{\rm{atoms}}>9$.

Considering the MAE for $E_{\rm SI}$ confirms these general findings, see
Figure LABEL:sifig:natomselementvsmae. For C-atoms the MAE for all reference
databases decreases for Set1 except for $n_{\rm atoms}=8$ for ANI-1E and
ANI-1, and increases moderately for Set2 with increasing value of $n_{\rm
atoms}$ for QM9 and ANI-1x and more steeply for PC9, ANI-1E and ANI-1 (see
Figures LABEL:sifig:natomselementvsmaeA and B). For nitrogen and oxygen atoms
satisfactory performance is only found for PC9, ANI-1x and ANI-1E (Set1) and
for QM9 and ANI-1x (Set2), see Figures LABEL:sifig:natomselementvsmaeC to F.
Molecules with a small number of atoms of a given element have a reduced
number of different chemical environments (see Figure
LABEL:sifig:NatomselementperDatabases). This makes it more difficult to
predict $\Delta E_{\rm Tauto}$ if that chemical environment is present in the
target data set (Tautobase) but not sampled in the reference sets.
Consequently, larger errors are observed for molecules with few atoms of a
given element.

###  4.2 Structural composition of the chemical databases

The structural diversity of the databases can also be quantified in terms of
the bond types that are covered. It can be assumed that the NN model learns
that specific composition, and consequently, if the database used for training
a NN model covers a large range of bond lengths, better results are expected.
In the following, bond length distributions in the reference databases of
equilibrium structures (PC9, QM9, and ANI-1E) are compared with the
distributions contained in Tautobase. Figures LABEL:sifig:dist_bonds_C_9 to
LABEL:sifig:dist_bonds_N_9p show that the reference and target distributions
have a different coverage of bond lengths. The general finding is that for
Set1 the overlap between reference and target distributions is better than for
Set2.

Figure LABEL:sifig:dist_bonds_C_9 shows that C-C single bonds between C(sp3)
atoms are well covered for the three reference databases compared with
Tautobase. The C(sp2)-C(sp2) double bonds are covered differently for the
reference datasets: QM9 has the fewest examples of this type of bond, whereas
ANI-1E shows the best coverage. Such bonds are important for large molecules
($n_{\rm{atoms}}>9$) because of the presence of aromatic rings (Figure
LABEL:sifig:dist_bonds_C_9p). Double C(sp2)-C(sp2) bonds close to hetero atoms
are poorly covered by all reference datasets. Those bonds are crucial because
they are the main origin of tautomerization rearrangement.

C(sp2)-N double bonds (Figure LABEL:sifig:dist_bonds_C_9) are abundantly
present in the Tautobase. However, the coverage of the reference datasets of
that type of bond is heterogeneous; ANI-1E shows the best coverage followed by
PC9 and QM9. On the other hand, C(sp2)-N bonds close to a heteroatom, more
prevalent in larger molecules, are better covered by QM9 than PC9 whereas
C(sp)-N bonds are well covered by all three databases. Carbon-Oxygen bonds for
carbonyl groups are more predominant in Set1 and are well covered for the
reference databases. Bonds for enols, esters and others are important for the
Tautobase; PC9 covers such C-O bonds sufficiently but it is poorly sampled for
QM9. Lastly, while C-O bonds of the type of alcohols and dialkyl ethers are
most sampled for the reference databases they are least important for the
Tautobase.

Figure 7: Values of the Kullback-Leibler divergence for different types of
bonds present in the reference databases (QM9,PC9 and ANI-1E) compared with
tautobase.

A quantitative measure for the overlap of two distributions is the KL
divergence $D(p||q)$, see Equation 2. If the two distributions are identical,
$D(p||q)=0$. On the other hand, if the reference database $p(x)$ (here QM9,
PC9, ANI-1E) contains more information than the target set $q(x)$ (Tautobase),
$D(p||q)>0$, and if specific information is missing, $D(p||q)<0$. Hence, the
cases for which $D(p||q)<0$ are of particular relevance if improvements of the
reference databases are sought for better capturing $\Delta E_{\rm Tauto}$.

The KL divergence analysis indicates that the coverage of the reference sets
is heterogeneous, see Tables LABEL:sitab:KL_9, LABEL:sitab:KL_9p and Figure 7.
There are several types of C-C bonds that are insufficiently covered, such as
C(sp2)-C(sp2) single and double bonds, or C(ar)-C(ar) bonds. Also, certain
types of C-N bonds would require more data as the bonds involving C(sp3) and
C(sp2) with different types of nitrogens. Coverage of C-O bonds by the
reference databases displays a bias toward alcohols, ethers and esters.
Finally, N-N are the types of bond that show a more diverse coverage between
databases, with some cases for which QM9 has a good coverage (N(3)-N(2) and
N(2)-N(2)(aromatic)) but a poor coverage for N(3)-N(3). Interestingly, there
are cases for which QM9 has a good coverage, whereas ANI-1E and PC9 are
deficient. Figure 7 shows that none of the reference databases covers all of
the predominant types of bonds present in the tautobase.

Next, the MAE for a specific number of a particular type of bond (e.g. C-C,
C-O, or C-N) was determined for single isomer energies, see Figure 8. The
results in Figure 8A show that for C-C bonds and Set1 the error for PC9
(orange) and ANI-1x (black) is constant and well below 1 kcal/mol. On the
other hand, for QM9 (red) the error oscillates without following a clear
trend. ANI-1E (light blue) and ANI-1 (dark blue) behave similarly to one
another with a smaller MAE for ANI-1E than the one for ANI-1.

For C-O bonds, the MAE of the prediction of $E_{\rm SI}$ slowly increases for
PC9 but remains well below 1 kcal/mol, whereas for QM9 it starts at above 5
kcal/mol and decreases to below 1 kcal/mol but always remaining above that for
PC9, see Figure 8C. The error for the database of the ANI family is largely
constant over the number of bonds. For ANI-1 the MAE oscillates between 1
kcal/mol and 2 kcal/mol, whereas for ANI-1E and ANI-1x the MAE is well below 1
kcal/mol, except for zero C-O bonds and ANI-1E. Considering C-N bonds (Figure
8E) it is found that their maximum number is larger than that for C-O bonds.
The magnitude for the MAE for this bond type is at least a factor of three
larger than that for the C-C and C-O bonds, respectively. Again, PC9 and
ANI-1x perform best, followed by ANI-1E (except for molecules with only one
C-N bond). The MAE for QM9 slowly decreases whereas that for ANI-1 is constant
at below 2 kcal/mol up to five C-N bonds after which it sharply increases.

Regarding C-O and C-N bonds, it is clear that the good coverage of ANI-1E
helps to reach small MAE when the number of bonds increases (Figure
LABEL:sifig:dist_bonds_C_9). These results show that PC9 has a good overall
performance because there is an adequate coverage of different chemical bond
types whereas QM9 and ANI-1E have biases toward some types of bonds (Table
LABEL:sitab:KL_9 and LABEL:sitab:KL_9p). It should be stressed though that
such an analysis excludes the fact that the same type of bond can behave
differently given different chemical environments.

For Set2 the increase of the error with decreasing number of samples is more
apparent. As discussed before, the MAE observed for larger molecules
($n_{\rm{atoms}}>9$) grows proportional with the number of bonds (See Figure 8
B, D and F). In this regard, ANI-1 and ANI-1E are the databases with largest
growth rate, followed by PC9 and finally QM9 and ANI-1x. The low MAE for
$E_{\rm SI}$ by number of bonds obtained with ANI-1x is a consequence of the
addition of an adequate number of non-equilibrium structures. This suggests
that a lack of chemical diversity can be partially compensated by including
non-equilibrium structures in a database.

Figure 8: Mean Absolute Error (MAE) of the prediction of the energy of single
isomers by number of bonds involving carbon atoms. The number of molecules for
different numbers of bonds is shown as a histogram. Panels A and B correspond
to the results of C-C bonds, panels C and D show the results for C-O bonds and
panels E and F show the results for C-N bonds. Left panel for Set1 and right
panel for Set2.

###  4.3 Initial Geometry

In the previous sections, the energy of the molecules was computed using the
equilibrium geometry of the tautomeric pairs computed at the level of theory
of the various reference databases used for training the NNs. However, in
practice, it would be of interest to sidestep the computationally rather
expensive optimization of the structures in the reference dataset (here
Tautobase) at the density functional or even higher level of quantum chemical
theory. For this, using empirical force fields is a possibility. A recent
study for 3271 small organic molecules ($n_{\rm atoms}<50$), similar to those
contained in Tautobase, found typical RMSDs of 0.25 Å to 1 Å between optimized
structures at the B3LYP/6-31G* level of theory compared with those optimized
by nine different force fields, including GAFF, MMFF94, OPLS and others.73
Considering this, it is interesting to assess the performance of the NN-based
models on FF-optimized geometries.

For this analysis, the geometries for molecules from Tautobase were optimized
as described in section 2.3 with the MMFF94 force field and then used to
evaluate the tautomerization energies using the five trained NNs. Table
LABEL:sitab:mae_MMFF94 shows that the MAE for the tautomerization and single
isomer energy increases for all evaluated models when the geometry used to
evaluate the energies differs from geometries optimized with the respective ab
initio method. In all cases the MAE for $\Delta E_{\rm Tauto}$ increases by a
factor between 1.5 and 3 compared with the error obtained using the optimized
geometries at the quantum chemical level of theory used to train the NN (see
Table 2). A similar effect is observed for $E_{\rm SI}$. It is noticeable that
this geometry effect is less pronounced for databases which contain non-
equilibrium structures: ANI-1 shows the smallest increase of the MAE for
$\Delta E_{\rm Tauto}$ compared with results from using optimized geometries
at the appropriate level of theory.

Normalized distributions for $\Delta E_{\rm Tauto}$ using MMFF94 geometries
for Set1 (Figure LABEL:sifig:error_dist_MMFF94C) indicate that the datasets
which only cover chemical space (QM9, PC9 and ANI-1E) perform similarly with
the highest values of the outliers close to 15 kcal/mol. Conversely, the other
two databases (ANI-1 and ANI-1x) have a more compact distribution with the
maximum values for outliers around 10 kcal/mol. Figure
LABEL:sifig:error_dist_MMFF94D shows that the most challenging case for
predicting tautomerization energies is for Set2 with geometries generated with
MMFF94. However, the performance is similar for all the datasets evaluated
with outliers larger than 20 kcal/mol, except for ANI-1x with a maximum of 10
kcal/mol. These results demonstrate that the initial geometry passed to the NN
model is essential for obtaining meaningful results. Scoring a model trained
on minimum energy structures computed at a given level of quantum chemical
theory can not be done using optimized structures from an empirical force
field (or from structures at a sufficiently different level of quantum
chemical theory).

To confirm this finding, a third set of molecules was evaluated using
geometries generated by six popular force fields as described in Section 2. In
the SAMPL2 challenge,61 the RMSD to DFT optimized geometries of several
molecules were evaluated with respect to geometries obtained from various
force fields. The energy predictions for the molecules on the SAMPL2 challenge
(Figure LABEL:sifig:RMSD_scatter) show no correlation between the geometry and
the energy predicted by the NN. There are molecules with a small RMSD (e.g.
$\leq 0.1$ Å) which display a significant error ($>5$ kcal/mol) in predicting
the energy by the NN models and vice versa. A possible explanation is that the
change in geometry can be compensated by more extensive sampling of chemical
space by the reference databases.

###  4.4 Visualization of chemical space

To understand the influence of the different databases studied on the
performance of the models, it is of interest to analyze the coverage of
‘chemical space’. Firstly, molecules in the databases were deconstructed and
their constituent amons (unique chemical fragments) were enumerated. PC9
contained the largest number of unique amons (8424), followed by QM9 (3929)
and, finally, the ANI family of databases (1663). There is significant overlap
of common amons between the datasets (Figure 9A, which suggests that they
cover similar regions of chemical space. Regarding the overlap of the test set
(Tautobase) with the databases tested, PC9 is the one which covers the most
amons by number in the reference set, followed by ANI-1E and, QM9 (Figure 9B).

Figure 9: A) Venn diagram showing the overlap of amons between the QM949,
PC950 and ANI51 family data sets. B) Overlap of amons up to length five
between three popular quantum chemistry datasets (QM9, PC9 and ANI-1E) and the
‘Tautobase’46, a collection of experimentally observed tautomers. Molecules
containing other atoms than hydrogen, carbon, nitrogen, and oxygen were
filtered from all datasets.

It is of interest to quantify the overlap of amons between the reference and
target sets as done in figure 9. This provides a measure for the coverage of
chemical space common to both sets and can be related to prediction errors for
$\Delta E_{\rm Tauto}$. Comparing this overlap with the results obtained in 2,
the database that covers the largest number of amons gives the smallest MAE
for $\Delta E_{\rm Tauto}$. This is more evident if the MAEs for $\Delta
E_{\rm Tauto}$ for Set1 and Set2 are analyzed separately. For Set1, the
results are consistent with the classification by number of amons (Figure 2 A,
C and E). However, for Set2 the results of QM9 are better than those obtained
with ANI-1E which contradicts the correlation with the number of amons in the
databases (Figure 2 B, D and F). This can probably be explained by the fact
that QM9 contains more large molecules which should help predicting the energy
of single molecules accurately. This effect is also observed looking at the
MAE for $E_{\rm SI}$ in Table 2 for which QM9 has a smaller MAE than PC9.

A more detailed analysis is possible by considering if the amons of the
isomers in the tautobase are present (or not) in the training databases. For
this, a set ‘seen amons’ (all constituent amons included in the reference
database) and ‘unseen amons’ (molecules for which one or several amons were
missing from the reference set) was defined. The error distributions for both
sets were determined and are reported in Figure 10. Perhaps unsurprisingly,
the ‘seen amons’ had a larger probability of obtaining a small error compared
to the ‘unseen amons’. Interestingly, PC9, which provides the broadest
sampling of chemical space as quantified by the number of amons in the
database, showed a similar probability error distribution for Set1 and Set2.
The errors for ‘unseen amons’ using the NN trained on QM9, a significantly
smaller dataset, shows a larger and more right-skewed distribution of errors.
One possible explanation may be that a better exploration of chemical space
helps when predicting energies for molecules containing chemistry outside that
covered by the database.

Regarding the ANI family of databases, the ANI-1 and ANI-1E results have
similar error distributions. However, ANI-1x shows a smaller mean error for
both the seen and unseen sets. These results are another indication that a
random sampling of conformational space does not help improve the NN model
predictions. On the contrary, it makes it worse than when only equilibrium
structures are considered. Another notable finding is that ANI-1x shows
similar performance for molecules with seen and unseen amons. This can be
explained given the good sampling of chemical space, which is the same as for
ANI-1E, but combined with a broad exploration of conformational space by a
refinement from ANI-1 using active learning 57.

Figure 10: The relationship between chemical space and model reliability.
Panels A to E: The mean average absolute error for molecules from Tautobase
with one or more amons outside the training set is larger than if all amons
are present, a trend observed for all databases. Panel F: Projection of the
Tautobase using the TMAP algorithm identifies ‘branches’ of chemical space
that are poorly predicted by the neural network models. The error displayed
here is for the QM9 database. TMAPs for all databases used for training are
available as interactive plots which can be viewed in a web browser, which can
be obtained in the supporting information.

Rational detection of systematic deficiencies in quantum chemical databases is
challenging because of the high dimensionality of chemical space. For this
reason, methods to visualise chemical space in a digestible way are highly
desirable. A recent development is the TreeMAP (TMAP) algorithm, which allows
for an interpretable, low dimensional representation of the test set’s
chemical space74. The TMAP algorithm constructs a weighted graph that is
efficient for compact representations of high dimensional data. The necessary
weights are based on the Jaccard distance, which measures the dissimilarity
between the fingerprint of two structures. This graph is then pruned to the
minimum spanning tree, a fully connected, acyclic sub-graph containing all
nodes of the parent graph, and retaining only essential edges which minimize
the weights. This organizes the compounds contained in a database into a tree,
putting related structures on nearby branches. From this, groups of related
moieties can be identified which are potentially detrimental to predicting the
quantity in question, here $\Delta E_{\rm Tauto}$.

For instance, coloring the nodes of Tautobase TMAP by error of tautomerization
energy (Figure 10F) reveals that structures with azoles containing N-N and N-O
bonds correlate with large errors. Interestingly, KL-divergences for these
types of bond distances suggested that they were underrepresented in the
reference sets, see Figure 7. The moieties corresponding to large errors
change based on the different databases used to train the NN model, see Figure
LABEL:sifig:TMAP_SI. Interactive plots are available in the supporting
information.

##  5 Discussion and Conclusions

The prominence of ML has raised concerns regarding the ’interpretability’ of
the models conceived75, 24. This awareness also increases for complex models
because a rational relationship between initial data used for training and
resulting prediction becomes less transparent. Therefore, it is important to
develop quantifiable and intuitive tests for how ML models “work” and how
trustworthy predictions by them are76. One recently proposed procedure is
“post-hoc” interpretation for which the practitioner analyzes a trained model
with the aim to understand what the model has learned from the data without
changing the underlying model23, 21. In this work, post-hoc interpretability
techniques were used to investigate the effect of different features of the
database on predicting a chemical property (tautomerization energy). The
selected features are considered important for the construction of robust
quantum chemical databases for ML. In the present case this implied the
analysis of individual features of several databases to quantify how these
modify the prediction of a chemical property on an unseen set of examples
using statistics and visualization techniques. With sufficient information
from such an analysis it is expected that it will be possible to identify
which features of the training databases are essential for good performance on
a given task, making a rational design/enhancement of databases for training
ML models for a given task possible.

The present work aimed at quantifying and analyzing the suitability of NNs
trained on five different reference data sets (QM9, PC9, ANI-1E, ANI-1, and
ANI-1x) to predict the tautomerization energies of molecules contained in
Tautobase. It was found that depending on which characteristics are
considered, the predicted MAEs can behave very differently and can, in part,
be related to geometrical and/or chemical properties encoded in the databases.
Such analyses attempt to digress from “black box” applications of ML methods
and move towards “interpretable ML”. Hence, one of the questions is “what
features need to be present and covered in a training database for application
to a concrete chemical question”. In the present case the databases to choose
from were QM9, PC9, ANI-1E, ANI-1, and ANI-1x’ and the application was
computation of the gas phase tautomerization energy.

The results indicate that the exploration of chemical space is essential for
meaningful results. The coverage of chemical space can be quantified by the
chemical diversity expressed as the number of amons on the database (see
Figure 9). The energy prediction improves when the overlap of the number of
amons in the training set and the tautobase increases. If the number of amons
in the chemical database does not cover all the amons on the target set,
addition of non-equilibrium structures to the training set can improve the
results. That addition needs to be done following a rational strategy because
an arbitrary addition causes deficiencies in the model, leading to poor
results as seen for the results for ANI-1E and ANI-1.

Another determinant property is the number of heavy atoms in molecules covered
in the database. Not surprisingly, better results are obtained for the range
covered by the database. Outside that range, the energy prediction quality
decreases with the number of atoms for most databases. One of the training
databases (ANI-1x) shows good results because the non-equilibrium structures
help in predicting the energies. It was observed that the different chemical
environments need to be thoroughly sampled because functional diversity is key
to assure good results.

The structural composition of the data sets used for training the NNs (QM9,
PC9 and ANI-1E) and the data set to which the trained models were applied to
(Tautobase) can be compared through the Kullback-Leibler divergence. The
overlap between these distributions already provides an indication how
suitable a particular reference data set will be for application to the target
task. In other words: the KL divergence can be used for the rational design of
databases for NN models. It will be of interest to extend this to angles and
dihedrals for a comprehensive exploration of the structural overlap. The
geometry of the molecule to evaluate the trained NN is essential for good
performance. Using the TMAP algorithm it was possible to identify regions of
chemical space that are poorly covered by the trained models. This is an
excellent aid because it becomes easier to improve sampling of a specific
region of chemical space to obtain better results. The characteristics of the
databases analyzed in the present work can be used as a rational basis to
determine whether a database is suitable for a specific prediction task.

In conclusion, the present work demonstrates that ML-trained models on five
different reference databases and applied to one specific task
(tautomerization energy) perform with a MAE ranging from 1.7 kcal/mol to 4.6
kcal/mol. The best performing reference database (ANI-1x with 5 M structures)
performs on average by 1 kcal/mol better than PC9 which contains about two
orders of magnitude fewer reference structures ($\approx 85$ K). On the other
hand, PC9 is chemically more diverse by a factor of 5 (as judged from the
number of amons) compared with the ANI family of databases. This indicates
that lack in chemical diversity can be compensated for by increased number of
non-equilibrium structures. However, the scaling of these two properties is
very different. Together with quantitative descriptors, such as the KL
divergence, the present results and analyses suggest that a rational approach
to database generation for specific tasks may be possible.

## Data Availability Statement

The PhysNet codes are available at Github
(https://github.com/MMunibas/PhysNet) . The database for ANI-1E
(10.5281/zenodo.4680953) and the geometries used for the tautobase
(10.5281/zenodo.4680972) are available at Zenodo.

## Acknowledgment

The authors acknowledge financial support from the Swiss National Science
Foundation (NCCR-MUST and Grant No. 200021-7117810) and the University of
Basel. OTU acknowledges funding from the Swiss National Science Foundation
(Grant No. P2BSP2_188147).

## References

  * Butler et al. 2018 Butler, K. T.; Davies, D. W.; Cartwright, H.; Isayev, O.; Walsh, A. Machine learning for molecular and materials science. _Nature_ 2018, _559_ , 547–555
  * Unke et al. 2021 In Press Unke, O. T.; Chmiela, S.; Sauceda, H. E.; Gastegger, M.; Poltavsky, I.; Schütt, K. T.; Tkatchenko, A.; Müller, K.-R. Machine learning force fields. _Chem. Rev._ 2021 In Press,
  * von Lilienfeld and Burke 2020 von Lilienfeld, O. A.; Burke, K. Retrospective on a decade of machine learning for chemical discovery. _Nat. Comm._ 2020, _11_ , 1–4
  * Noé et al. 2020 Noé, F.; Tkatchenko, A.; Müller, K.-R.; Clementi, C. Machine learning for molecular simulation. _Annu. Rev. Phys. Chem._ 2020, _71_ , 361–390
  * Agrawal and Choudhary 2016 Agrawal, A.; Choudhary, A. Perspective: Materials informatics and big data: Realization of the “fourth paradigm” of science in materials science. _Apl Mater._ 2016, _4_ , 053208
  * Aspuru-Guzik et al. 2018 Aspuru-Guzik, A.; Lindh, R.; Reiher, M. The matter simulation (r) evolution. _ACS Cent. Sci._ 2018, _4_ , 144–152
  * Behler and Parrinello 2007 Behler, J.; Parrinello, M. Generalized neural-network representation of high-dimensional potential-energy surfaces. _Phys. Rev. Lett._ 2007, _98_ , 146401
  * Bartók et al. 2010 Bartók, A. P.; Payne, M. C.; Kondor, R.; Csányi, G. Gaussian approximation potentials: The accuracy of quantum mechanics, without the electrons. _Phys. Rev. Lett._ 2010, _104_ , 136403
  * Rupp et al. 2012 Rupp, M.; Tkatchenko, A.; Müller, K.-R.; Von Lilienfeld, O. A. Fast and accurate modeling of molecular atomization energies with machine learning. _Phys. Rev. Lett._ 2012, _108_ , 058301
  * Montavon et al. 2013 Montavon, G.; Rupp, M.; Gobre, V.; Vazquez-Mayagoitia, A.; Hansen, K.; Tkatchenko, A.; Müller, K.-R.; Von Lilienfeld, O. A. Machine learning of molecular electronic properties in chemical compound space. _New J. Phys._ 2013, _15_ , 095003
  * Faber et al. 2017 Faber, F. A.; Hutchison, L.; Huang, B.; Gilmer, J.; Schoenholz, S. S.; Dahl, G. E.; Vinyals, O.; Kearnes, S.; Riley, P. F.; Von Lilienfeld, O. A. Prediction errors of molecular machine learning models lower than hybrid DFT error. _J. Chem. Theory Comput._ 2017, _13_ , 5255–5264
  * Schütt et al. 2018 Schütt, K. T.; Sauceda, H. E.; Kindermans, P.-J.; Tkatchenko, A.; Müller, K.-R. SchNet–A deep learning architecture for molecules and materials. _J. Chem. Phys._ 2018, _148_ , 241722
  * Unke and Meuwly 2018 Unke, O. T.; Meuwly, M. A reactive, scalable, and transferable model for molecular energies from a neural network approach based on local information. _J. Chem. Phys._ 2018, _148_ , 241708
  * Wilkins et al. 2019 Wilkins, D. M.; Grisafi, A.; Yang, Y.; Lao, K. U.; DiStasio, R. A.; Ceriotti, M. Accurate molecular polarizabilities with coupled cluster theory and machine learning. _Proc. Natl. Acad. Sci. USA_ 2019, _116_ , 3401–3406
  * Veit et al. 2020 Veit, M.; Wilkins, D. M.; Yang, Y.; DiStasio Jr, R. A.; Ceriotti, M. Predicting molecular dipole moments by combining atomic partial charges and atomic dipoles. _J. Chem. Phys._ 2020, _153_ , 024113
  * Unke et al. 2020 Unke, O. T.; Koner, D.; Patra, S.; Käser, S.; Meuwly, M. High-dimensional potential energy surfaces for molecular simulations: from empiricism to machine learning. _Mach. Learn.: Sci. Technol._ 2020, _1_ , 013001
  * Ko et al. 2021 Ko, T. W.; Finkler, J. A.; Goedecker, S.; Behler, J. A fourth-generation high-dimensional neural network potential with accurate electrostatics including non-local charge transfer. _Nat. Comm._ 2021, _12_ , 1–11
  * Von Lilienfeld 2018 Von Lilienfeld, O. A. Quantum machine learning in chemical compound space. _Angew. Chem. Int. Ed._ 2018, _57_ , 4164–4169
  * Heinen et al. 2020 Heinen, S.; Schwilk, M.; von Rudorff, G. F.; von Lilienfeld, O. A. Machine learning the computational cost of quantum chemistry. _Mach. Learn.: Sci. Technol._ 2020, _1_ , 025002
  * Käser et al. 2020 Käser, S.; Koner, D.; Christensen, A. S.; von Lilienfeld, O. A.; Meuwly, M. Machine Learning Models of Vibrating H2CO: Comparing Reproducing Kernels, FCHL, and PhysNet. _J. Phys. Chem. A_ 2020, _124_ , 8853–8865
  * Du et al. 2019 Du, M.; Liu, N.; Hu, X. Techniques for interpretable machine learning. _Commun. ACM_ 2019, _63_ , 68–77
  * Samek and Müller 2019 Samek, W.; Müller, K.-R. _Explainable AI: interpreting, explaining and visualizing deep learning_ ; Springer, 2019; pp 5–22
  * Murdoch et al. 2019 Murdoch, W. J.; Singh, C.; Kumbier, K.; Abbasi-Asl, R.; Yu, B. Definitions, methods, and applications in interpretable machine learning. _Proc. Natl. Acad. Sci. USA_ 2019, _116_ , 22071–22080
  * Dybowski 2020 Dybowski, R. Interpretable machine learning as a tool for scientific discovery in chemistry. _New J Chem_ 2020, _44_ , 20914–20920
  * Wilkinson and McNaught 1997 Wilkinson, A.; McNaught, A. IUPAC Compendium of Chemical Terminology,(the" Gold Book"). _International Union of Pure and Applied Chemistry: Zürich, Switzerland_ 1997,
  * Raczyńska et al. 2005 Raczyńska, E. D.; Kosińska, W.; Ośmiałowski, B.; Gawinecki, R. Tautomeric equilibria in relation to pi-electron delocalization. _Chem Rev_ 2005, _105_ , 3561–3612
  * Martin 2009 Martin, Y. C. Let’s not forget tautomers. _J. Comput. Aided Mol. Des._ 2009, _23_ , 693–704
  * Watson and Crick 1953 Watson, J. D.; Crick, F. H. Molecular structure of nucleic acids: a structure for deoxyribose nucleic acid. _Nature_ 1953, _171_ , 737–738
  * Shukla and Leszczynski 2013 Shukla, M. K.; Leszczynski, J. Tautomerism in nucleic acid bases and base pairs: a brief overview. _Wires Comput. Mol Sci._ 2013, _3_ , 637–649
  * Singh et al. 2015 Singh, V.; Fedeles, B. I.; Essigmann, J. M. Role of tautomerism in RNA biochemistry. _RNA_ 2015, _21_ , 1–13
  * Käser et al. 2020 Käser, S.; Unke, O. T.; Meuwly, M. Isomerization and decomposition reactions of acetaldehyde relevant to atmospheric processes from dynamics simulations on neural network-based potential energy surfaces. _J. Chem. Phys._ 2020, _152_ , 214304
  * Sitzmann et al. 2010 Sitzmann, M.; Ihlenfeldt, W.-D.; Nicklaus, M. C. Tautomerism in large databases. _J. Comput. Aided Mol. Des._ 2010, _24_ , 521–551
  * Greenwood et al. 2010 Greenwood, J. R.; Calkins, D.; Sullivan, A. P.; Shelley, J. C. Towards the comprehensive, rapid, and accurate prediction of the favorable tautomeric states of drug-like molecules in aqueous solution. _J. Comput. Aided Mol. Des._ 2010, _24_ , 591–604
  * Pospisil et al. 2003 Pospisil, P.; Ballmer, P.; Scapozza, L.; Folkers, G. Tautomerism in computer-aided drug design. _J. Recept. Signal Transduct._ 2003, _23_ , 361–371
  * Taylor et al. 2014 Taylor, P. J.; van der Zwan, G.; Antonov, L. Tautomerism: Introduction, history, and recent developments in experimental and theoretical methods. _Tautomerism: methods and theories_ 2014, 1–24
  * Fogarasi 2010 Fogarasi, G. Studies on tautomerism: benchmark quantum chemical calculations on formamide and formamidine. _J. Mol. Struct_ 2010, _978_ , 257–262
  * Käser et al. 2020 Käser, S.; Unke, O. T.; Meuwly, M. Reactive dynamics and spectroscopy of hydrogen transfer from neural network-based reactive potential energy surfaces. _New J. Phys._ 2020, _22_ , 055002
  * Baughcum et al. 1981 Baughcum, S. L.; Duerst, R. W.; Rowe, W. F.; Smith, Z.; Wilson, E. B. Microwave Spectroscopic Study Of Malonaldehyde (3-Hydroxy-2-Prop enal) .2. Structure, Dipole-Moment, and Tunneling. _J. Am. Chem. Soc._ 1981, _103_ , 6296–6303
  * Firth et al. 1991 Firth, D. W.; Beyer, K.; Dvorak, M. A.; Reeve, S. W.; Grushow, A.; Leopold, K. R. Tunable far infrared spectroscopy of malonaldehyde. _J. Chem. Phys._ 1991, _94_ , 1812–1819
  * Limbach et al. 1987 Limbach, H.-H.; Wehrle, B.; Zimmermann, H.; Kendrick, R. D.; Yannoni, C. S. Kinetic 15N-CPMAS-NMR Study of a Double Proton Transfer in a Crystalline Malonaldehyde Diimine Derivative. _Angew. Chem. Int. Ed._ 1987, _26_ , 247–248
  * Wang et al. 2008 Wang, Y.; Braams, B. J.; Bowman, J. M.; Carter, S.; Tew, D. P. Full-dimensional quantum calculations of ground-state tunneling splitting of malonaldehyde using an accurate ab initio potential energy surface. _J. Chem. Phys._ 2008, _128_ , 224314
  * Schröder et al. 2011 Schröder, M.; Gatti, F.; Meyer, H.-D. Theoretical studies of the tunneling splitting of malonaldehyde using the multiconfiguration time-dependent Hartree approach. _J. Chem. Phys._ 2011, _134_ , 234307
  * Yang and Meuwly 2010 Yang, Y.; Meuwly, M. A generalized reactive force field for nonlinear hydrogen bonds: Hydrogen dynamics and transfer in malonaldehyde. _J. Chem. Phys._ 2010, _133_ , 064503
  * Karandashev et al. 2017 Karandashev, K.; Xu, Z.-H.; Meuwly, M.; Vaníček, J.; Richardson, J. O. Kinetic isotope effects and how to describe them. _Struct. Dyn._ 2017, _4_ , 061501
  * Unke and Meuwly 2019 Unke, O. T.; Meuwly, M. PhysNet: a neural network for predicting energies, forces, dipole moments and partial charges. _J. Chem. Theory Comput._ 2019, _15_ , 3678–3693
  * Wahl and Sander 2020 Wahl, O.; Sander, T. Tautobase: An Open Tautomer Database. _J. Chem. Inf. Model._ 2020, _60_ , 1085–1089
  * Dhaked et al. 2020 Dhaked, D. K.; Guasch, L.; Nicklaus, M. C. Tautomer Database: A Comprehensive Resource for Tautomerism Analyses. _J. Chem. Inf. Model._ 2020, _60_ , 1090–1100
  * Huang and von Lilienfeld 2020 Huang, B.; von Lilienfeld, O. A. Quantum machine learning using atom-in-molecule-based fragments selected on the fly. _Nat. Chem._ 2020, 1–7
  * Ramakrishnan et al. 2014 Ramakrishnan, R.; Dral, P. O.; Rupp, M.; Von Lilienfeld, O. A. Quantum chemistry structures and properties of 134 kilo molecules. _Sci. Data_ 2014, _1_ , 140022
  * Glavatskikh et al. 2019 Glavatskikh, M.; Leguy, J.; Hunault, G.; Cauchy, T.; Da Mota, B. Dataset’s chemical diversity limits the generalizability of machine learning predictions. _J. Cheminf._ 2019, _11_ , 69
  * Smith et al. 2017 Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1, A data set of 20 million calculated off-equilibrium conformations for organic molecules. _Sci. Data_ 2017, _4_ , 3192–3203
  * Smith et al. 2020 Smith, J. S.; Zubatyuk, R.; Nebgen, B.; Lubbers, N.; Barros, K.; Roitberg, A. E.; Isayev, O.; Tretiak, S. The ANI-1ccx and ANI-1x data sets, coupled-cluster and density functional theory properties for molecules. _Sci. Data_ 2020, _7_ , 1–10
  * Ruddigkeit et al. 2012 Ruddigkeit, L.; Van Deursen, R.; Blum, L. C.; Reymond, J.-L. Enumeration of 166 billion organic small molecules in the chemical universe database GDB-17. _J. Chem. Inf. Model._ 2012, _52_ , 2864–2875
  * Nakata and Shimazaki 2017 Nakata, M.; Shimazaki, T. PubChemQC project: a large-scale first-principles electronic structure database for data-driven chemistry. _J. Chem. Inf. Model._ 2017, _57_ , 1300–1308
  * Fink et al. 2005 Fink, T.; Bruggesser, H.; Reymond, J.-L. Virtual exploration of the small-molecule chemical universe below 160 daltons. _Angew. Chem. Int. Ed._ 2005, _44_ , 1504–1508
  * Fink and Reymond 2007 Fink, T.; Reymond, J.-L. 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. _J. Chem. Inf. Model._ 2007, _47_ , 342–353
  * Smith et al. 2018 Smith, J. S.; Nebgen, B.; Lubbers, N.; Isayev, O.; Roitberg, A. E. Less is more: Sampling chemical space with active learning. _J. Chem. Phys._ 2018, _148_ , 241733
  * Smith et al. 2017 Smith, J. S.; Isayev, O.; Roitberg, A. E. ANI-1: an extensible neural network potential with DFT accuracy at force field computational cost. _Chem. Sci._ 2017, _8_ , 3192–3203
  * O’Boyle et al. 2011 O’Boyle, N. M.; Banck, M.; James, C. A.; Morley, C.; Vandermeersch, T.; Hutchison, G. R. Open Babel: An open chemical toolbox. _J. Cheminf._ 2011, _3_ , 33
  * Halgren 1996 Halgren, T. A. Merck molecular force field. I. Basis, form, scope, parameterization, and performance of MMFF94. _J. Comput. Chem._ 1996, _17_ , 490–519
  * Geballe et al. 2010 Geballe, M. T.; Skillman, A. G.; Nicholls, A.; Guthrie, J. P.; Taylor, P. J. The SAMPL2 blind prediction challenge: introduction and overview. _J. Comput. Aided Mol. Des._ 2010, _24_ , 259–279
  * Foloppe and MacKerell 2000 Foloppe, N.; MacKerell, A. D., Jr All-atom empirical force field for nucleic acids: I. Parameter optimization based on small molecule and condensed phase macromolecular target data. _J. Comput. Chem._ 2000, _21_ , 86–104
  * Wang et al. 2004 Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and testing of a general amber force field. _J. Comput. Chem._ 2004, _25_ , 1157–1174
  * Jorgensen and Tirado-Rives 1988 Jorgensen, W. L.; Tirado-Rives, J. The OPLS force field for proteins. Energy minimizations for crystals of cyclic peptides and crambin. _J. Am. Chem. Soc._ 1988, _110_ , 1657–1666
  * Rappé et al. 1992 Rappé, A. K.; Casewit, C. J.; Colwell, K.; Goddard III, W. A.; Skiff, W. M. UFF, a full periodic table force field for molecular mechanics and molecular dynamics simulations. _J. Am. Chem. Soc._ 1992, _114_ , 10024–10035
  * Schmid et al. 2011 Schmid, N.; Eichenberger, A. P.; Choutko, A.; Riniker, S.; Winger, M.; Mark, A. E.; van Gunsteren, W. F. Definition and testing of the GROMOS force-field versions 54A7 and 54B7. _Eur. Biophys. J._ 2011, _40_ , 843–856
  * Hassinen and Peräkylä 2001 Hassinen, T.; Peräkylä, M. New energy terms for reduced protein models implemented in an off-lattice force field. _J. Comput. Chem._ 2001, _22_ , 1229–1242
  * Stewart 2013 Stewart, J. J. Optimization of parameters for semiempirical methods VI: more modifications to the NDDO approximations and re-optimization of parameters. _J. Mol. Model._ 2013, _19_ , 1–32
  * Stewart 2016 Stewart, J. J. MOPAC2016. _Stewart Computational Chemistry: Colorado Springs, CO, USA_ 2016,
  * Frisch et al. 2009 Frisch, M. J.; Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. et al. Gaussian 09, Revision D. 01, Gaussian. _Inc.: Wallingford, CT_ 2009,
  * Diwekar and David 2015 Diwekar, U.; David, A. _BONUS Algorithm for Large Scale Stochastic Nonlinear Programming Problems_ ; Springer, 2015; pp 27–34
  * Cover and Thomas 2006 Cover, T. M.; Thomas, J. A. _Elements of Information Theory_ ; Wiley Series in Telecommunications and Signal Processing; John Wiley & Sons, 2006
  * Lim et al. 2020 Lim, V. T.; Hahn, D. F.; Tresadern, G.; Bayly, C. I.; Mobley, D. L. Benchmark assessment of molecular geometries and energies from small molecule force fields. _F1000Res._ 2020, _9_ , 1390
  * Probst and Reymond 2020 Probst, D.; Reymond, J.-L. Visualization of very large high-dimensional data sets as minimum spanning trees. _J. Cheminf._ 2020, _12_ , 12
  * Lipton 2018 Lipton, Z. C. The Mythos of Model Interpretability: In machine learning, the concept of interpretability is both important and slippery. _ACM Queue_ 2018, _16_ , 31–57
  * Schütt et al. 2019 Schütt, K. T.; Gastegger, M.; Tkatchenko, A.; Müller, K.-R. _Explainable AI: Interpreting, Explaining and Visualizing Deep Learning_ ; Springer, 2019; pp 311–330

Table of Content