#  Ab initio machine learning of phase space averages

Jan Weinreich  University of Vienna, Faculty of Physics, Kolingasse 14-16,
AT-1090 Wien, Austria  University of Vienna, Vienna Doctoral School in
Physics, Boltzmanngasse 5, 1090 Vienna, Austria  Dominik Lemm  University of
Vienna, Faculty of Physics, Kolingasse 14-16, AT-1090 Wien, Austria
University of Vienna, Vienna Doctoral School in Physics, Boltzmanngasse 5,
1090 Vienna, Austria  Guido Falk von Rudorff  University of Vienna, Faculty of
Physics, Kolingasse 14-16, AT-1090 Wien, Austria  O. Anatole von Lilienfeld
<EMAIL_ADDRESS>Machine Learning Group, Technische
Universität Berlin, 10587 Berlin, Germany  Berlin Institute for the
Foundations of Learning and Data, 10587 Berlin, Germany.  Institute of
Physical Chemistry and National Center for Computational Design and Discovery
of Novel Materials (MARVEL), Department of Chemistry, University of Basel,
Klingelbergstrasse 80, CH-4056 Basel, Switzerland

###### Abstract

Equilibrium structures determine material properties and biochemical
functions. We propose to machine learn phase-space averages, conventionally
obtained by ab initio or force-field based molecular dynamics (MD) or Monte
Carlo simulations. In analogy to ab initio molecular dynamics (AIMD), our ab
initio machine learning (AIML) model does not require bond topologies and
therefore enables a general machine learning pathway to ensemble properties
throughout chemical compound space. We demonstrate AIML for predicting
Boltzmann averaged structures after training on hundreds of MD trajectories.
AIML output is subsequently used to train machine learning models of free
energies of solvation using experimental data, and reaching competitive
prediction errors (MAE $\sim$ 0.8 kcal/mol) for out-of-sample molecules–
within milli-seconds. As such, AIML effectively bypasses the need for MD or
MC-based phase space sampling, enabling exploration campaigns throughout CCS
at a much-accelerated pace. We contextualize our findings by comparison to
state-of-the-art methods resulting in a Pareto plot for the free energy of
solvation predictions in terms of accuracy and time.

##  I Introduction

Figure 1:  Example of a chemical compound space (CCS) for $N_{\rm p}=N_{\rm
e}=32$ protons (and electrons) with a hierarchy of composition, constitution,
and conformation. Each level corresponds to distinct temperature-regimes and
is described by specific quantum chemical methods1; 2; 3; 4; 5; 6; 7; 8; 9;
10. Ab initio machine learning (AIML) can act on all three levels and does not
require fixed constitutions but allows general ensemble predictions
$\mathcal{A}$ using ML-based averaged structures $\mathbfcal{R}$ and free
energy machine learning11 (FML).

Structure determines function – a hallmark paradigm in the atomistic sciences,
ranging from biologists studying protein functions based on x-ray structures
to organic chemists discussing reaction mechanisms based on NMR
measurements12. The connection between structure $\mathbfcal{R}$ and a
compound’s function is given through statistical mechanics averages
$\mathcal{A}$ over the ensemble $\mathcal{E}$ of Boltzmann weighted
configurations,

$\displaystyle\mathcal{E}\rightarrow\mathbfcal{R}\rightarrow\mathcal{A}~{}.$
(1)

The function of a molecule depends on the biological or chemical context e.g.
solubilities or binding affinities – all of which can be expressed as phase
space averages $\mathcal{A}$. Understanding this relation is of fundamental
importance as the temperature-dependent balance of configurations dictates the
biological function of proteins and their macroscopic behavior (think of egg-
white). Unfortunately, to quantitatively predict thermal averages which
minimize the free energy imposes major computational challenges due to the
necessity of sampling phase space. Furthermore, since experimental efforts to
obtain a compounds’ structure $\mathbfcal{R}$ are cumbersome several
computational routes have been introduced.

However, covering molecular structures poses a monumental challenge. We also
highlight the bigger picture of chemical compound space13 (CCS) and its
hierarchical structure given by composition, constitution, and conformation.
The inherent curse of the dimensions of CCS means that even considering all
possible molecules of a single fixed composition quickly results in a
combinatorial explosion as illustrated in Fig. 1. Thus most approaches follow
a divide-and–and-conquer strategy addressing the combinatorial problem of CCS
at individual levels. Still, the numerical complexity of studying such
relationships using molecular dynamics14 (MD) or Monte Carlo15; 16 (MC) is
overwhelming and to this day most methods with accurate, yet rapid predictions
suffer from the curse of conformer sampling. For instance, atomistic
simulations study statistical mechanics (SM) ensembles through molecular
dynamics
$\mathcal{E}\xrightarrow{\text{MD}}\mathbfcal{R}\xrightarrow{\text{SM}}\mathcal{A}$
and are deeply intertwined with insights into biological functions. Ab initio
molecular dynamics (AIMD) simulations not only allow studying molecules but
also chemical reactions2; 17; 18; 19. However, they are much more costly than
force fields20; 5; 6; 21; 1 due to having to solve approximate quantum
mechanical equations at every time step. To this account hybrid set-ups
$\mathcal{E}\xrightarrow{\text{MD}}\mathbfcal{R}\xrightarrow{\text{ML}}\mathcal{A}$
using both atomistic simulation and machine learning (ML) have been introduced
uniting quantum mechanical equations with surrogate learning on the fly
potentials22; 23; 24. This helps mitigate some of the ab initio costs but may
still require extensive MD sampling. These challenges have driven
technological advancements of dedicated computer hardware, e.g. of the
supercomputer Anton25, specifically designed to accelerate MD simulations.
Conversely, large MD codes5; 26 have been rewritten in CUDA27 just so that
they can run on GPUs. Decentralized global computing network initiatives such
as Folding$@$home28; 29 also predominantly run MD. MD also routinely consumes
major fractions of resources and energy costs of high-performance computing
centers, as recently seen for the Gordon Bell award to Car and co-workers for
running MD on 100 M atoms111https://www.acm.org/media-
center/2020/november/gordon-bell-prize-2020.

Figure 2:  Centered histograms and standard deviation $\sigma_{d}$ of all
distances between ten non-hydrogen atoms of aspirin extracted from an ab
initio molecular dynamics (AIMD) trajectory31 (a). Conventional AIMD and ab
initio ML (AIML) map the ensemble $\mathcal{E}$ to averaged structure
$\mathbfcal{R}$ to the statistical mechanics (SM) average $\mathcal{A}$ (b).
Bold labeled index pairs (i,j) define a bond topology (c). Sketch of two
principal components $\text{PC}1,\text{PC}2$ of ML-based representations32;
33; 34; 35 on a fictitious free energy surface (d) corresponding to conformers
represented by a disconnectivity graph36. Instead of MD average
representations11 (FML) $\mathbf{X}$ we propose AIML predicted representations
of averaged conformers $\mathbf{x}(\mathbfcal{R})$.

To address the length and timescale problem of conformational space across
CCS, we have recently introduced Free energy Machine Learning (FML) which
relies on the averaged structure $\mathbfcal{R}$ as input to predict ensemble
averages such as the free energy 11,
$\mathbfcal{R}\xrightarrow{\text{FML}}\mathcal{A}$. Averaging the structure is
necessary because ensemble properties inherently depend on multiple
configurations and as such using a single geometry per molecule introduces
ambiguities to the ML model. Here, we propose to replace the preceding step,
i.e. the generation of the averaged input through extensive molecular dynamics
runs for any query compound by an ab initio machine learning (AIML) model,
$\mathcal{E}\xrightarrow{\text{AIML}}\mathbfcal{R}$.

AIML makes use of the Graph-To-Structure37 (G2S) method to predict three-
dimensional structures though chemical compound space. Training data and
labels, however, are fundamentally different from G2S: Instead of using a
single optimized conformer per molecule, AIML training data consists of
averages over complete MD trajectories, enabling the prediction of
thermodynamically averaged conformers. This accounts for the fundamentally
important difference between a Boltzmann average and a single atomic
configuration for ensemble property predictions, as previously discussed in
free energy machine learning11 (FML).

In particular, we replace ensemble sampling with machine learning of ensemble
averages $\mathcal{A}$ of an equilibrium property $a(\mathbf{r},\mathbf{p})$
by,

$\displaystyle\mathbfcal{A}=\;\frac{1}{Z}\int
a(\mathbf{r},\mathbf{p})~{}e^{-\beta
E(\mathbf{r},\mathbf{p})}~{}\text{d}\mathbf{r}~{}\text{d}\mathbf{p}\approx\mathbfcal{A}^{\rm
ML}(\mathbfcal{R})~{},$ (2)

with $\beta$, $Z$, and $E$ being the Boltzmann-factor, the partition function,
and the total energy, respectively38. The approximate equality of Eq. 2 is
achieved by training an ML model $\mathbfcal{A}^{\text{ML}}$ that uses only
the averaged conformer $\mathbfcal{R}$ with values $\mathbfcal{A}$ that
include the rigorous integral over the ensemble.

AIML is a purely ML-based framework and properly accounts for the underlying
Boltzmann statistics and can subsequently be used to generate the appropriate
input for FML model-based predictions. The goal of our work is not to build an
ML model that perfectly reflects thermodynamic expectation values but to
construct a surrogate model that can predict these integrals with high
accuracy and speed. By including Boltzmann averaged conformers we make sure to
include a canonical mapping of the underlying ensemble of each molecule to the
ensemble property. By training a second FML model on experimental values we
ensure that our mapping $\mathbfcal{A}^{\text{ML}}(\mathbfcal{R})$ from the
averaged structure also includes contributions from the complete ensemble as
defined in the proper phase space integral $\mathbfcal{A}$. Our numerical
results, i.e. the systematic improvement of the models’ accuracy with training
set size, indicate that our assumption (Eq. 5) is sufficiently valid for the
dataset we studied.

After briefly introducing AIML in the following, we will demonstrate its
applicability for the prediction of aqueous free energies of solvation of out-
of-sample molecules without having to perform explicit MD simulations. For
training, however, extensive MD trajectories at corresponding temperatures are
necessary, as well as experimental measurements of solvation reference
energies. Lastly, we provide an overview of the efficiency of AIML in the
context of alternative state-of-the-art solvation methods. We find that AIML
offers respective speed-ups by four to seven orders of magnitude when compared
to classical or ab initio MD-based predictions of free energies of solvation.

##  II Methods

###  II.1 FML and G2S

For the previously published free energy machine learning11 (FML) approach
first all sampled geometries $\mathbf{r}$ had to be transformed to
representation vectors $\mathbf{x}(\mathbf{r})$ before finally the average
$\mathbf{X}$ over the representation vectors,

$\displaystyle\mathbf{X}=\int\mathbf{x}(\mathbf{r})~{}e^{-\beta
E(\mathbf{r})}~{}\text{d}\mathbf{r}~{},$ (3)

could be computed. The key advantage of AIML is that instead of explicitly
sampling conformer space, a single AIML evaluation is required to predict a
surrogate vector $\mathbf{x}(\mathbfcal{R})$ evaluated for the system average
$\mathbfcal{R}$ to replace $\mathbf{X}$ (s. Fig. 2d).

Graph-To-Structure (G2S) exploits implicit correlations among relaxed
structures in training data sets to infer interatomic distances for out-of-
sample compounds across chemical space. G2S effectively enables direct
reconstruction of three-dimensional coordinates, thereby allows circumventing
conventional energy optimization. G2S can reach an accuracy on par or better
than conventional structure generators. As query input, G2S requires only
bond-network and stoichiometry-based information. G2S learns the direct
mapping from a chemical graph to that structure that had been recorded in the
training data set. For the prediction of new structures, only molecular
connectivity is needed, which can be provided e.g. via SMILES39 or SELFIES40.
The G2S machines predict all pairwise distances. The full 3D geometry is then
reconstructed using DGSOL41 for heavy atoms and a Lebedev sphere optimization
scheme for hydrogen atoms.

###  II.2 Workflow

Figure 3:  To generate ab initio machine learning (AIML) training data,
conformer space is sampled to obtain Boltzmann weighted average distance
matrices for given molecular graphs. AIML predicts the averaged distance
matrix, which is then converted to a three-dimensional geometry
$\mathbfcal{R}$. Finally, the ensemble average $\mathcal{A}$ is predicted
using free energy machine learning11 (FML).

An essential ingredient for AIML was extending the previous Graph-To-
Structure37 (G2S) method via the introduction of an averaged geometry
$\mathbfcal{R}$. This enables a computationally efficient ML-based map between
ensemble and free energies and addresses the conformer sampling bottleneck.
AIML provides an effective representation for the conformer ensemble by
mapping all degrees of freedom to a single averaged structure to approximate
the ensemble-averaged structure with corresponding AIML representation vector
$\mathbf{x}$ (cf. Fig. 2). A schematic overview of the steps for AIML training
and prediction is given in Fig. 3. As we will discuss in the following two
steps are needed to combine both methods: i) use Boltzmann weighted distance
matrices $\mathbf{D}$ for training ii) use AIML average conformer predictions
as input for an ensemble average model (s. Fig. 3 right). The required
training steps for the AIML structure prediction use Boltzmann-averaged
intramolecular distance matrices. More specifically, as a first step, the
molecule is transformed to a graph-based representation37 $\mathbf{g}$ with
the average distance matrix $\mathbf{D}$ as training labels resulting in ML
models for heavy atom pairs as well as heavy and hydrogen atoms. Before
entering the next step of free energy prediction, AIML is trained with the
maximal number of molecules ($N=512$) to construct the average training and
test set conformers. This process is repeated for the complete dataset with
consistent training test splits between AIML structure and free energies
prediction. Next, a machine (s. sec. II.3) for learning free energies is
trained using AIML predicted average geometries. Finally, AIML can be used for
out-of-sample predictions (s. Fig. 3). Based on the molecular graph, the
Boltzmann weighted distance matrix $\mathbf{D}$ is predicted. Next, a
distance-geometry solver41 (DGSOL) is used to convert the distance matrix to
three-dimensional coordinates. The predicted average conformer then serves as
a link between the graph and corresponding three-dimensional geometry (s. Fig.
2). Subsequently, the average predicted conformer $\mathbfcal{R}$ is
transformed to a single Bag-of-Bonds42 (BoB) representation vector
$\mathbf{x}(\mathbfcal{R})$ and used to predict the free energy.

###  II.3 Machine learning

ML, in particular of solutions to the Schrodinger equation i.e. quantum
machine learning43; 44; 45; 46; 47 (QML), allows navigating chemical compound
space (CCS) with high efficiency and precision. ML has become a popular avenue
to material science with applications to atomization energies48; 49, crystal
formation energies50, carbenes51, excited states52; 53, oxidation states54,
nuclear magnetic resonance spectra55, reaction barriers56; 57, magnetic
systems58 and charge transfer59 or molecular fragments60; 61. Many recent
approaches are based on individual levels of CCS, i.e. composition62; 63,
constitution64; 65; 66; 67 (s. Fig. 1). Moreover ensemble properties such as
the free energy of solvation11; 68; 69; 70; 71; 72; 73; 70; 74, melting
points24; 75; 76, magnetic anisotropy tensors77, phases of water78; 79; 80
have previously been addressed with ML.

Our ML approach is based on Kernel-Ridge Regression81 (KRR) a supervised
learning method that allows approximating arbitrary functional relationships
between input data given as molecular representations $\mathbf{x}$ and
properties $\mathcal{A}(\mathbf{x})$. Using KRR $\mathbf{x}$ is mapped into a
high dimensional feature space rendering the regression problem linear. A
remarkable result of KRR81; 82; 83 is that the mapping does not need to be
carried out explicitly, instead, the distances between representations
$\mathbf{x}_{i}$ and $\mathbf{x}_{j}$ are computed e.g. using Gaussian kernel
functions,

$\displaystyle
k(\mathbf{x}_{i},\mathbf{x}_{j})=\exp{\left(-\frac{||\mathbf{x}_{i}-\mathbf{x}_{j}||^{2}_{2}}{2\sigma^{2}}\right)}~{},$
(4)

that measure the similarity between two compounds $i$ and $j$ resulting in the
kernel matrix $\mathbf{K}$ where $\sigma$ is the kernel-width hyperparameter.
We use KRR to predict the vector of all interatomic distances $\mathbf{D}_{q}$
of a query compound $q$ and a distance geometry solver41 (DSGOL) for
subsequent reconstruction of the three-dimensional geometry.

The vector $\mathbf{D}_{q}$ contains all distances of atoms in the query
molecule. The average conformer prediction
$\mathbfcal{R}(\mathbf{x}_{q}=\mathbf{g}_{q})$ of a query molecule $q$
represented by a graph37 $\mathbf{g}_{q}$ is given by,

$\displaystyle\mathbfcal{R}(\mathbf{g}_{q})=\text{DGSOL}(\mathbf{D}(\mathbf{g}_{q}))~{},$
(5)

where the distance prediction is as follows,

$\displaystyle\mathbf{D}_{q}=\mathbf{K}\cdot\boldsymbol{\alpha}~{}.$ (6)

Here, the kernel matrix $\mathbf{K}$ is evaluated between query and training
compounds $\mathbf{g}_{i}$ and $\mathbf{g}_{q}$ with regression coefficients
$\boldsymbol{\alpha}$. The optimal regression coefficients
$\boldsymbol{\alpha}$ are obtained by solving a set of equations,

$\displaystyle\boldsymbol{\alpha}=(\mathbf{K}+\lambda\cdot\mathds{I})^{-1}\mathbf{D}~{},$
(7)

where the vector $\mathbf{D}$ contains all distances between atoms for each of
the training molecules and $\lambda$ is a regularization parameter. The
ensemble property prediction is given by:

$\displaystyle\mathbfcal{A}_{q}=\widetilde{\mathbf{K}}\cdot\widetilde{\boldsymbol{\alpha}}~{}.$
(8)

As before, for training $\widetilde{\mathbf{K}}$ is evaluated between all
training compounds now using molecular representation vectors
$\mathbf{x}(\mathbfcal{R})$:

$\displaystyle\widetilde{\boldsymbol{\alpha}}=(\tilde{\mathbf{K}}+\lambda\cdot\mathds{I})^{-1}\mathbfcal{A}~{},$
(9)

where the vector $\mathbfcal{A}$ contains the values of the ensemble property
in the training set. Learning curves quantify the model prediction error,
often measured as mean absolute error (MAE), against the number of training
samples $N$ and are key to understand the efficiency of ML models. It is
generally found81 that they are linear on a log-log scale,

$\displaystyle\log{\left(\frac{\text{MAE}}{\text{unit}}\right)}\approx
I-S\cdot\log{(N)}~{},$ (10)

where $I$ is the initial error and $S$ is the slope indicating model
improvement given more training data.

##  III Results

###  III.1 Concept of equilibrium structure prediction

Within AIML, we view the averaged structure as an ensemble property,
$\mathbfcal{R}=\mathcal{A}$, representing the connection to the aforementioned
overarching theme of structure determining function12. In addition, the
thermal equilibrium structure is relevant for NMR spectroscopy which accounts
for protein flexibility by resulting in time-averaged structures equivalent to
$\mathbfcal{R}$ due to the ergodic theorem38. Purely ML-based implementation
bypasses routinely encountered sampling issues – de facto replacing MD
simulations with predicted averages as ensemble fingerprints for subsequent
property prediction.

We first use kernel-ridge regression 81 (KRR) to predict the symmetric matrix
of averaged interatomic distances of a query compound $q$, i.e.
$\mathbf{D}_{q}\approx\mathbf{K}\cdot\boldsymbol{\alpha}$ which contains all
inferred averaged distances of atoms in the query molecule. $\mathbf{K}$ and
$\boldsymbol{\alpha}$ correspond to the kernel matrix and training weights
obtained for a training set consisting of MD trajectories and averaged
interatomic distances as labels. In analogy to Graph-To-Structure 37 (G2S), we
subsequently rely on the distance geometry solver41 (DSGOL) for reconstruction
of the three-dimensional structure $\mathbfcal{R}$, as well as on graph-based
representations, $\mathbf{g}_{q}$. To exemplify the AIML approach, consider
the ab initio molecular dynamics (AIMD) trajectory published in Ref. 31 of the
aspirin molecule at 300 K, resulting atomic distance histograms in Fig. 2a).
In order to establish a graph-based representation to replace ab initio MD
with ab initio ML (AIML) as illustrated in Fig. 2b, it is necessary to assign
bonds. This is straightforward using the distance histogram as covalently
bonded atoms will not move far from each other (see Fig. 2c)) – a concept that
can be generalized via coarse-graining21; 84; 85; 86. While AIML requires a
suggested molecular graph it is not restricted to a single fixed bond topology
but allows adapting the molecular graph depending on the relevant degrees of
freedom depending on temperature or the environment of the molecule. Thus AIML
can include the formation and breaking of bonds (cfg. Fig. 1) corresponding to
adding or erasing a one in the bond topology matrix $\mathbf{g}_{q}$. (cfg.
Fig. 2). Note that the predicted pairwise distance matrix does not only
account for nearest neighbor effects but the complete many-body description
since it includes all cross combinations of atomic distances.

Using the graph as the representation for constructing kernels, the AIML model
then learns the center of each off-diagonal element in the distance histogram
as a label. Next, the AIML predicted distances are used to reconstruct the
average conformer. Subsequently the ML representation vector
$\mathbf{x}(\mathbfcal{R})$ is computed. Therefore, AIML allows exchanging the
order of average evaluation compared to the previous FML11 approach
(illustrated Fig. 2) resulting in a dramatic reduction of computational costs.

AIML proposes a different paradigm by connecting all hierarchies of CCS with
ensemble properties into a single ML-based framework. Because of its
generality, AIML does not require a priori information about bonds but only
averaged atomic distances. In analogy to the hierarchy of CCS, AIML has
several special cases with fundamental physical interpretations (s. Fig. 1):
At temperatures higher than most bond energies AIML operates on atomic
clusters corresponding to a topology matrix that (mostly) contains zeros i.e.
the composition. At moderate temperatures, bonds exist but may occasionally
break corresponding to adding or removing a zero in the topology matrix. In
this case, multiple molecular graphs can be extracted and AIML predicts the
constitution averaged structure.

###  III.2 Application: Ensemble to structure to property

Figure 4:  Learning curve of the Root-Mean-Square Deviation (RMSD), a measure
of structural distance, of ab initio machine learning (AIML) three-dimensional
structure predictions as a function of the number of training structures
$N_{s}$ using GAFF25; 6 force field and density functional theory (DFT) with
Becke-Perdew87; 88 (BP) functional for training structure sampling (a). Two
predicted average conformers $\mathbfcal{R}$ show the improvement of structure
prediction along the learning curve (at $N_{s}=32,128,512$). The mean absolute
error (MAE) of predicted free energies $\mathcal{F}$ of the FreeSolv89
database as a function of $N$ for free energy machine learning (FML) (b) using
MD sampling11 versus AIML (no sampling) and the Bag-of-Bonds42 (BoB)
representation. The AIML $N_{s}=32$ model was trained with only 32 structures
for conformer prediction.

In this section, we demonstrate the usefulness of AIML for the problem of
accurate predictions of free energies of solvation. In particular, we focus on
experimental free energies of solvation of 642 charge-neutral small to medium-
sized bio-organic molecules, as encoded in the FreeSolv database 89. Solvation
free energies90; 5; 6; 91; 92; 93; 94; 95; 96; 97; 98; 99; 100 are of
fundamental importance for chemistry, and the FreeSolv database has become a
popular benchmark for performance testing of novel models68; 69; 70; 71; 72;
73; 70; 74; 24; 75.

To use AIML to predict averaged structures and FML to predict free energies of
solvation for out-of-sample molecules, we first trained AIML models as
described above using molecular graphs as input and as labels for the averaged
distances.

These were obtained from extended force-field based MD runs and density
functional theory (DFT) for conformer sampling with Def2TZVPD-FINE101; 102;
103; 104 basis set and Becke-Perdew87; 88 (BP) functional in the gas-phase (s.
sec. IV.1 for details). Note that we neglect the effect of water on the phase
space of the solute (e.g. through hydrogen bonds). We believe that this aspect
warrants further in-depth investigation within subsequent studies in the
future.

Depending on the temperature (cfg. Fig. 1), some degrees of freedom do not get
averaged out by the phase space integral (Eq. 2), restricting the domain to
certain local basins of the total free energy. These remaining degrees of
freedom can be identified in the way described in Fig 2a). We demonstrate the
idea of AIML for ambient temperatures in the gas phase for which conventional
molecular graph topologies as in biochemistry hold – without any loss of
generality. Within such a regime, we can safely assume that any distance
histogram matrix would be consistent with one topology which represents the
coarse-grained back-bone that is not averaged out by the phase space integral.
Of course, this approach could also be applied to any other set of conditions
presuming that there is some way to easily infer valid topologies as a
function of conditions (like temperature and composition). The latter is a
separate problem that goes beyond the scope of this work.

As numerical evidence of the functionality of the AIML idea, we present in
Fig. 3a) prediction errors of the three-dimensional Boltzmann averaged
structures $\mathbfcal{R}$ as a function of training set size (aka learning
curves81; 105; 106).

Numerical results shown in Fig. 4 a) indicate a systematic linear improvement
on a log-log scale107 as a function of the size of the training set $N_{s}$,
that is, the number of averaged training structures. Note that $N$ stands for
the number of training points for free energy values and that we have trained
and evaluated two different machines after each other, for structure and free
energy prediction respectively. For the maximal training set size considered
($N_{s}=512$ molecules, 80% of FreeSolv), the average root-mean-square
deviation108; 109; 110 (RMSD), a measure of structural distance between
structures, has decayed to only 0.80 Å for predicted FF and 0.82 Å for DFT
averaged conformers, the slope of the learning curve, however, indicates that
learning has not yet been saturated. We find structure prediction of large
molecules a particularly hard problem i.e. on average the RMSDs increase with
the size of the structure. This was also observed for a random subset of the
GDB17 molecular database 111 (s. SI. Fig. 7b for the scatter plot of molecular
size vs. RMSD indicating a rough correlation).

Note that the corresponding learning curve predicting optimal (not Boltzmann-
averaged) distances is less steep for the molecules in FreeSolv, and exhibits
a higher offset (s. SI Fig. 4). This could be due to the fact that learning
thermal averages is less ambiguous for AIML than learning potential energy
minima as it is in the case of G2S. This also holds for the heavy atom-
hydrogen distances (s. SI Tab. 1 and Fig. 1 in SI). From a different point of
view, by computing the Boltzmann average over distances conformer flexibility
is effectively integrated out, thus simplifying learning compared to the
optimized structures.

Next, we demonstrate the learning efficiency (s. Fig. 4b) of AIML for free
energy prediction based on previously predicted structures $\mathbfcal{R}$.
The main disadvantage of the preceding free energy machine learning11 (FML)
model was that it requires explicit conformer sampling for free energy
prediction. The novel advantage of AIML is that no sampling is required.
Instead, after prediction of the average conformer $\mathbfcal{R}$ the free
energy prediction is based on the ML representation vector
$\mathbf{x}(\mathbfcal{R})$ where $\mathbf{x}$ is the ML-based representation
Bag-of-Bonds42 (BoB). Encouragingly, FML (requiring explicit MD sampling) and
AIML exhibit similar learning curves, achieving mean absolute errors (MAE) of
$0.68\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$ and
$0.82\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$, respectively, as
shown in Fig. 4b). This indicates that the performance of distance-based
representations in conjunction with AIML is fairly robust showing only
$0.14\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$ loss of accuracy
compared to running MD simulations. Using predicted averaged DFT conformers we
obtain roughly the same accuracy of
$0.84\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$ (s. SI. Fig. 8). This
is consistent with our previous assessment of out-of-sampling AIML conformer
RMSDs using FF and DFT which also resulted in comparable errors for both
methods. To illustrate the importance of structure prediction for subsequent
property prediction, we have added a learning curve in Fig. 4b) using an AIML
model with only $N=32$ average conformers for training the structure
prediction. The mentioned model shows a much smaller learning rate. Generally,
we find that better average structure prediction will lead to improved
subsequent FML models (discussed in more detail in SI in Fig. 5).

We find AIML to perform consistently better when the training distances result
from Boltzmann sampling instead of using optimized structures - similar to
what we also noted above for the structure prediction (see SI Fig. 9). Even
more surprising is the observation that training free energy prediction with
AIML predicted structures resulted in slightly better models than training
with the ground-truth averaged conformers resulting from MD. This indicates
that AIML effectively smoothens conformer space by isolating the most
important degrees of freedom, thus facilitating structure-based regression of
thermal averages.

###  III.3 Assessment of efficiency

Figure 5:  Comparison of solvation methods (if multiple weighted conformers
used red else blue) in terms of mean absolute error (MAE) and order of
magnitude of per molecule prediction time $t$ for FreeSolv89 database. Pareto
front (dotted) is formed by methods with best accuracy and cost per prediction
trade-off. Central processing unit (CPU) compute time varies depending on
hardware, code, etc., and was estimated if not available. All references and
MAE of free energies are listed in the SI Tab. 2.

To gain a more comprehensive idea of AIML’s value to the field, we have
assessed the cost accuracy trade-off. Testing AIML on the FreeSolv89 database
we have measured average prediction times of
$41\text{\,}\mathrm{ms}$/molecule. These prediction times are dominated by the
structure reconstruction task ($40\text{\,}\mathrm{ms}$), while only
$1\text{\,}\mathrm{ms}$ is required to yield the free energy estimate (on a
single-core AMD EPYC 7402P compute chip). For comparison, the corresponding
prediction based on a classical force-field MD simulation protocol would have
consumed three to four orders of magnitude more time, not to mention the costs
associated with quantum chemistry based ab initio MD. To gain a comprehensive
overview of the field, we have performed solvation free energy calculations
for all of the FreeSolv molecules using the following methods (all listed free
energies available, s. sec. V, MAEs listed in SI. Tab. II):

  1. 1.

Solvation model based on density94 (SMD) at M06-2X112/Def2-TZVPP101; 102; 103;
104 (timing for SMD implemented in Gaussian113 not to be published)

  2. 2.

COSMO-RS-B1 and COSMO-RS-B2 referring to COSMO-RS95; 96; 114; 115; 116 with
Def2TZVP101; 102; 103; 104 and Def2TZVPD-FINE101; 102; 103; 104 basis sets
respectively and Becke-Perdew87; 88 (BP) functional in the gas-phase

  3. 3.

Reaction mechanism generator group (RMG) solvation117

  4. 4.

ddCOSMO118 results obtained with PySCF119 and PBE-0120; 121/Def2TZVPD101; 102;
103; 104

  5. 5.

Generalized Born122; 123 (GBSA) model, results obtained using AMBER5; 6

  6. 6.

Free energy machine learning11 (FML) with explicit conformer sampling on
FreeSolv database

To complete the picture, we also included literature values for FreeSolv
concerning the methods ARROW-PIMD820, Thermodynamic integration (TI) with
GAFF25; 6 extracted from the FreeSolv89 database and reference interaction
site model124; 125; 100; 124; 126 (3D-RISM). The trade-off between cost and
accuracy, including an outline of the resulting Pareto front, is displayed in
Fig. 5.

We note that AIML adds to the convexity of the Pareto front, representing a
meaningful compromise: Although roughly twice as slow, it is slightly more
accurate at $0.82\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$ than the
Reaction Mechanism Generator127 (RMG) model (MAE of
$0.98\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$), but still four
orders of magnitude faster than the next best ab initio method COSMO-RS95; 96;
114; 115; 116. A list of all MAE is provided in SI Tab. 2. Thus, AIML is
positioned on the Pareto front of the available solvation methods located in a
sweet spot between speed and accuracy, providing the fastest predictions at
the given accuracy of about
$0.82\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$. Note that the AIML
learning curves have not yet saturated and that its accuracy will likely
further improve if more training samples are included (c.f. Fig. 4). Improving
the AIML model will hardly worsen the prediction time due to the linear
scaling of KRR predictions w.r.t. training set size (s. sec. II.3) and
therefore shift the Pareto front towards higher accuracy. Furthermore, it is
important to note that arbitrary accurate ab initio trajectories can be used
for training while the prediction time is independent of the level of theory.
We expect that recently published ML models tailored towards solvation such as
SoluteML128 may outperform the presented AIML models’ accuracy, but we note
that only a small training set of $N=512$ molecules was used and we expect MAE
to decay further with the training set size. A3D-PNAConv-FT129 combines the 2D
and 3D structure and transfer learning achieving an MAE of
$0.417\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$ for the FreeSolv
data set. Predictions require conformer sampling using a FF and the lowest-
energy conformer. In contrast, AIML does not require sampling and is in fact
replacing the functionality of a force-field or of an ab-initio calculation.

Note that we have also tried to combine RMG and AIML/FML via the $\Delta$-ML
approach130 where RMG is used as a baseline for AIML, but unfortunately, the
prediction errors did not improve (s. SI Fig. 8). Moreover, AIML performs
worse for large molecules with many conformers (s. SI Fig. 7b). Unfortunately,
combining a random sampled GDB17111 dataset with 10000 molecules with the
FreeSolv average conformers did not lead to improved structure predictions due
to the small overlap of the two data sets (s. SI Fig. 6). Specifically, the
small training set size and very high chemical diversity of the FreeSolv
database including the elements C, H, O, S, N, F, I, Br, P, Cl, and up to 24
non-hydrogen atoms per molecule limit the accuracy of structure prediction for
large compounds. Instead of adding random structures to improve structure
prediction for the FreeSolv database (s. SI Fig. 6) we could show that
sampling131 the local chemical space of the largest FreeSolv can help to
improve the models’ accuracy (s. SI Fig. 8). Alternatively, these problems may
be resolved by improved graph-based representations that include information
about local chemical substructures, leading to better structure and improved
free energy predictions.

Recent graph-based ML models132; 133 can achieve a competitive accuracy with
root mean squared errors (RMSEs) around
$1\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$. We have achieved a
similar RMSE of $1.35\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$ for a
training set size of $N=512$. We emphasize that the AIML approach is very
different: First, AIML uses three-dimensional conformations which can lead to
a much-improved accuracy (MAE of
$0.57\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$ for $N=490$) as we
have shown earlier11 and allows going beyond fixed graphs. Secondly, AIML also
predicts ensemble-based representation whereas SMILES-based ML use molecular
graphs as input. We found a direct comparison with the two previously
mentioned other ML methods132; 133 difficult because they either use a
different training-test split132 or neglect133 certain molecules of the
FreeSolv database89. The comparison shows that our method might have a
slightly higher initial offset due to having to learn the representation i.e.
the averaged conformer before predicting the free energy. Our learning curves
(s. Fig. 4) do not indicate saturation of the MAE with training set size $N$
and might still surpass graph-based models for large training set sizes
because AIML contains information about molecular conformations.

##  IV Conclusion

We have introduced ab initio machine learning (AIML) allowing for efficient
predictions of ensemble averages which systematically improve in accuracy as
training set sizes grow. To the best of our knowledge, for the first time,
AIML effectively bypasses the need for extensive MD or MC simulations to
directly infer Boltzmann averaged geometries. Unlike all other solvation
models (shown in Fig. 5) the AIML framework could easily be applied to other
ensemble properties, e.g. melting points, without much adaptation since no
manual pre-selection of features for molecular fingerprints is required. AIML
does not require any additional sampling for inferring ensemble averages of
new out-of-sample query molecules: Instead AIML accounts for multiple
Boltzmann weighted configurations implicitly through its training data. We
have exemplified AIML for estimating experimental solvation free energies, and
our numerical results amount to evidence showing that the conformer ensemble
can effectively be linked to a single averaged conformer that serves as a
canonical representative. AIML predictions are consistent with the previous
free energy machine learning11 (FML) approach without the need to run an MD
simulation for each prediction, reaching errors as low as
$0.82\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$ for 41 CPU-
ms/molecule prediction cost.

Further analysis has revealed that AIML does not yet work well with all
available molecular representations32. More specifically, we find that
representations, tailored toward atomization energies and including explicit
angular dependencies, such as FCHL19134; 135, yield less favorable AIML models
(s. SI Fig. 9). Conversely, it might be possible to further improve AIML by
tailoring and optimizing representations and architecture (e.g. using
locality, symmetry, neural networks).

The question of uniqueness is very fundamental for molecular
representations136; 137; 138; 139; 140; 141. For the present dataset, however,
this was not an issue as all averaged representations distinguished all data
items. It is possible, however, to imagine a scenario where this is not the
case: For two different ensembles, $\mathcal{E}$ and $\bar{\mathcal{E}}$ with
the same average conformer $\mathbfcal{R}=\bar{\mathbfcal{R}}$, but two
different averages $\mathcal{A}\neq\bar{\mathcal{A}}$, the presented AIML
model would make the same predictions and could not distinguish between these
two systems. However, this problem could be resolved by including higher-order
moments for the prediction of the representation, like including an AIML model
for the standard deviation of the representation. Future work will deal with
this question. Here our main focus was on molecules with uniquely defined bond
topologies. For future applications, AIML could be applied to cases where
assigning bond topologies is ambiguous or impossible such as transition
states142; 143 or molecules at very high temperatures. These are important
cases where AIML can be used, but graph-based ML models cannot be used.

In summary, in comparison to classical or ab initio MD-based predictions of
free energies of solvation, AIML offers respective speed-ups by four to seven
orders of magnitude. AIML achieves such speed-ups by effectively shifting
computational cost for the query prediction to the training set generation.
However, in light of the sheer scale of the chemical compound space available
for molecular queries, this trade-off might be useful.

###  IV.1 Conformer and Free Energy Data

ML based on a single geometry can lead to ambiguous predictions11 ensemble
property predictions because predictions can vary substantially depending on
the conformer. Many body representations134; 135; 144; 145; 146 rely on three-
dimensional geometries, which becomes even more relevant if the target
property depends on multiple relevant conformers. A solution to this issue is
to sample configuration space to obtain a conformer invariant ML
representations11. Sampling can be achieved by different strategies: MD
simulations5; 6, systematic conformer ensemble scans147; 148; 149 and
conformer generation methods either knowledge or force field (FF) based such
as ETKDG6150, Gen3D 151 and others152; 153; 154; 155; 156. A more expensive
method but accurate method is to obtain conformations using ab initio
approaches such as density functional theory (DFT)18; 19 or tight binding147;
148; 149 (TB). Despite these advantages, a common pitfall of these methods is
that sometimes extension to arbitrary chemistries is not straightforward. To
this end, ML-based methods37; 157; 158; 159 hold the promise of providing
faster and more general structure predictions. There are only very few ML
structure generation methods e.g., based on reinforcement learning160 or
stochastic normalizing flows161 that take energy weights of different
conformers into account. Here, to obtain a diverse set of conformers as the
AIML training set we have performed MD simulations in vacuum at an elevated
temperature of $T=$350\text{\,}\mathrm{K}$$ using OpenMM26 using a Langevin
integrator. GAFF25; 6 with a time-step of $\Delta t=$2\text{\,}\mathrm{fs}$$
was used with a total simulation time of $2\text{\,}\mathrm{ns}$. Partial
charges are computed with antechamber5; 6 at AM1-BCC162 level. MD samples are
selected with $2\text{\,}\mathrm{ps}$ time separation. To compare AIML with
COSMO-RS95; 96; 114; 115; 116 solvation method, we used the COSMO-RS workflow
based on ab initio DFT calculations with Turbomole163 and the Becke-Perdew
(BP)87; 88 functional as implemented in COSMOconf164 with two different basis
sets, Def2TZVP and Def2TZVPD-FINE101; 102; 103; 104 (for future reference
referred to as B1 and B2). Based on these results, free energies are extracted
using the COSMOtherm165 program. In addition, the reaction mechanism generator
group (RMG) based approach was used to compute free energies of solvation117
of the FreeSolv database via the leruli.com API166. The FreeSolv89 dataset
contains 642 charge neutral compounds and their experimental free energies of
solvation. The average unsigned error of the experimental values is
$0.57\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$ close to the level of
thermal energy fluctuations
($k_{B}\cdot$300\text{\,}\mathrm{K}$\approx$0.6\text{\,}\mathrm{kcal}\text{\,}{\mathrm{mol}}^{-1}$$).
All ML models use a maximal training set size of $80\%$ corresponding to
$N=512$ molecules. Hyperparameters are optimized with nested five-fold cross
validation.

##  V Data and Code availability

The AIML code and all free energies of solvation of the FreeSolv database (if
produced by the authors) are published in a freely available repository
https://doi.org/10.5281/zenodo.6401711. We gladly provide more data for
specific requests.

##  VI Acknowledgements

We acknowledge support from the European Research Council (ERC-CoG Grant QML).
This project has received funding from the European Union’s Horizon 2020
research and innovation program under Grant Agreement #772834. This research
was also supported by the NCCR MARVEL, a National Centre of Competence in
Research, funded by the Swiss National Science Foundation (grant number
182892).

##  VII Declaration of Conflicting Interests

DL, GFvR, and OAvL are shareholders of Leruli GmbH.

## References

  * van Duin  _et al._ (2001) A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard, The Journal of Physical Chemistry A 105, 9396 (2001).
  * Car and Parrinello (1985) R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
  * Oostenbrink  _et al._ (2004) C. Oostenbrink, A. Villa, A. E. Mark, and W. F. Van Gunsteren, Journal of Computational Chemistry 25, 1656 (2004), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.20090 .
  * Christen  _et al._ (2005) M. Christen, P. H. Hünenberger, D. Bakowies, R. Baron, R. Bürgi, D. P. Geerke, T. N. Heinz, M. A. Kastenholz, V. Kräutler, C. Oostenbrink, C. Peter, D. Trzesniak, and W. F. van Gunsteren, Journal of Computational Chemistry 26, 1719 (2005), https://onlinelibrary.wiley.com/doi/pdf/10.1002/jcc.20303 .
  * Wang  _et al._ (2006) J. Wang, W. Wang, P. A. Kollman, and D. A. Case, J. Mol. Graph. Model. 25, 247 (2006).
  * Wang  _et al._ (2004) J. Wang, R. M. Wolf, J. W. Caldwell, and P. A. Kollman, J. Comp. Chem. 25, 1157 (2004).
  * Vega  _et al._ (2009) C. Vega, J. L. F. Abascal, M. M. Conde, and J. L. Aragones, Faraday Discuss. 141, 251 (2009).
  * Jorgensen  _et al._ (1996) W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, Journal of the American Chemical Society 118, 11225 (1996).
  * Jorgensen  _et al._ (1983) W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, The Journal of Chemical Physics 79, 926 (1983).
  * Abrams  _et al._ (2006) J. B. Abrams, L. Rosso, and M. E. Tuckerman, The Journal of chemical physics 125, 074115 (2006).
  * Weinreich  _et al._ (2021) J. Weinreich, N. J. Browning, and O. A. von Lilienfeld, The Journal of Chemical Physics 154, 134113 (2021).
  * Voet and Voet (2006) D. Voet and J. G. Voet,  _Biochemistry_ (John Wiley and Sons, Inc., Hoboken, NJ, USA, 2006).
  * von Lilienfeld  _et al._ (2020) O. A. von Lilienfeld, K.-R. Müller, and A. Tkatchenko, Nature Reviews Chemistry 4, 347 (2020).
  * Frenkel and Smit (2001) D. Frenkel and B. Smit,  _Understanding Molecular Simulation_ , 2nd ed. (Academic Press, Inc., USA, 2001).
  * Rosenbluth (2003) M. N. Rosenbluth, AIP Conference Proceedings 690, 22 (2003).
  * Hastings (1970) W. K. Hastings, Biometrika 57, 97 (1970).
  * Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Kresse and Joubert (1999) G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Pereyaslavets  _et al._ (2022) L. Pereyaslavets, G. Kamath, O. Butin, A. Illarionov, M. Olevanov, I. Kurnikov, S. Sakipov, I. Leontyev, E. Voronina, T. Gannon, G. Nawrocki, M. Darkhovskiy, I. Ivahnenko, A. Kostikov, J. Scaranto, M. Kurnikova, S. Banik, H. Chan, M. Sternberg, and B. Fain, Nature Communications 13 (2022), 10.1038/s41467-022-28041-0.
  * Kmiecik  _et al._ (2016) S. Kmiecik, D. Gront, M. Kolinski, L. Wieteska, A. E. Dawid, and A. Kolinski, Chemical Reviews 116, 7898 (2016).
  * Csányi  _et al._ (2004) G. Csányi, T. Albaret, M. C. Payne, and A. De Vita, Phys. Rev. Lett. 93, 175503 (2004).
  * Li  _et al._ (2015) Z. Li, J. R. Kermode, and A. De Vita, Phys. Rev. Lett. 114, 096405 (2015).
  * Jinnouchi  _et al._ (2019) R. Jinnouchi, F. Karsai, and G. Kresse, Physical Review B 100, 014105 (2019).
  * Shaw  _et al._ (2008) D. E. Shaw, M. M. Deneroff, R. O. Dror, J. S. Kuskin, R. H. Larson, J. K. Salmon, C. Young, B. Batson, K. J. Bowers, J. C. Chao, M. P. Eastwood, J. Gagliardo, J. P. Grossman, C. R. Ho, D. J. Ierardi, I. Kolossváry, J. L. Klepeis, T. Layman, C. McLeavey, M. A. Moraes, R. Mueller, E. C. Priest, Y. Shan, J. Spengler, M. Theobald, B. Towles, and S. C. Wang, Commun. ACM 51, 91–97 (2008).
  * Eastman  _et al._ (2017) P. Eastman, J. Swails, J. D. Chodera, R. T. McGibbon, Y. Zhao, K. A. Beauchamp, L.-P. Wang, A. C. Simmonett, M. P. Harrigan, C. D. Stern, R. P. Wiewiora, B. R. Brooks, and V. S. Pande, PLOS comp. bio. 13, e1005659 (2017).
  * Vingelmann  _et al._ (2020) P. Vingelmann, F. Fitzek, and NVIDIA, NVIDIA https://developer. nvidia. com/cuda-toolkit (2020).
  * Zagrovic  _et al._ (2002) B. Zagrovic, C. D. Snow, M. R. Shirts, and V. S. Pande, Journal of Molecular Biology 323, 927 (2002).
  * Voelz  _et al._ (2010) V. A. Voelz, G. R. Bowman, K. Beauchamp, and V. S. Pande, Journal of the American Chemical Society 132, 1526 (2010).
  * Note (1) https://www.acm.org/media-center/2020/november/gordon-bell-prize-2020.
  * Chmiela  _et al._ (2017) S. Chmiela, A. Tkatchenko, H. E. Sauceda, I. Poltavsky, K. T. Schütt, and K.-R. Müller, Science Advances 3, e1603015 (2017).
  * Parsaeifard  _et al._ (2021a) B. Parsaeifard, D. S. De, A. S. Christensen, F. A. Faber, E. Kocer, S. De, J. Behler, O. A. von Lilienfeld, and S. Goedecker, Machine Learning: Science and Technology 2, 015018 (2021a).
  * Bartók  _et al._ (2013) A. P. Bartók, R. Kondor, and G. Csányi, Phys. Rev. B 87, 184115 (2013).
  * Nigam  _et al._ (2022) J. Nigam, G. Fraux, and M. Ceriotti, “Unified theory of atom-centered representations and graph convolutional machine-learning schemes,” (2022), arXiv:2202.01566 [stat.ML] .
  * Huang and von Lilienfeld (2016a) B. Huang and O. A. von Lilienfeld, The Journal of Chemical Physics 145, 161102 (2016a).
  * Wales (2012) D. Wales, Philosophical transactions. Series A, Mathematical, physical, and engineering sciences 370, 2877 (2012).
  * Lemm  _et al._ (2021a) D. Lemm, G. F. von Rudorff, and O. A. von Lilienfeld, Nature Communications 12, 1 (2021a).
  * Tuckerman (2001) M. Tuckerman,  _Statistical Mechanics: Theory And Molecular Simulation_ (Oxford Graduate Texts, 2001).
  * Weininger (1988) D. Weininger, Journal of Chemical Information and Computer Sciences 28, 31 (1988).
  * Krenn  _et al._ (2020) M. Krenn, F. Häse, A. Nigam, P. Friederich, and A. Aspuru-Guzik, Machine Learning: Science and Technology 1 (2020), 10.1088/2632-2153/aba947.
  * Moré and Wu (1999) J. J. Moré and Z. Wu, Journal of Global Optimization 15, 219 (1999).
  * Hansen  _et al._ (2015) K. Hansen, F. Biegler, R. Ramakrishnan, W. Pronobis, O. A. von Lilienfeld, K.-R. Müller, and A. Tkatchenko, The Journal of Physical Chemistry Letters 6, 2326 (2015).
  * Huang and von Lilienfeld (2020a) B. Huang and O. A. von Lilienfeld, “Ab initio machine learning in chemical compound space,” (2020a), arXiv:2012.07502 [physics.chem-ph] .
  * Ceriotti  _et al._ (2021) M. Ceriotti, C. Clementi, and O. Anatole von Lilienfeld, The Journal of Chemical Physics 154, 160401 (2021).
  * Westermayr  _et al._ (2021) J. Westermayr, M. Gastegger, K. T. Schütt, and R. J. Maurer, The Journal of Chemical Physics 154, 230903 (2021).
  * Tkatchenko (2020) A. Tkatchenko, Nature Communications 11, 1 (2020).
  * Keith  _et al._ (2021) J. A. Keith, V. Vassilev-Galindo, B. Cheng, S. Chmiela, M. Gastegger, K.-R. Müller, and A. Tkatchenko, Chemical Reviews 121, 9816 (2021).
  * Rupp  _et al._ (2012a) M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, Phys. Rev. Lett. 108, 058301 (2012a).
  * Faber  _et al._ (2017) F. A. Faber, L. Hutchison, B. Huang, J. Gilmer, S. S. Schoenholz, G. E. Dahl, O. Vinyals, S. Kearnes, P. F. Riley, and O. A. von Lilienfeld, J. Chem. Th. Comp. 13, 5255 (2017).
  * Faber  _et al._ (2016) F. A. Faber, A. Lindmaa, O. A. von Lilienfeld, and R. Armiento, Phys. Rev. Lett. 117, 135502 (2016).
  * Schwilk  _et al._ (2020) M. Schwilk, D. N. Tahchieva, and O. A. von Lilienfeld, “Large yet bounded: Spin gap ranges in carbenes,” (2020), arXiv:2004.10600 [physics.chem-ph] .
  * Westermayr and Marquetand (2021) J. Westermayr and P. Marquetand, Chemical Reviews 121, 9873 (2021).
  * Dral and Barbatti (2021) P. O. Dral and M. Barbatti, Nature Reviews Chemistry 5, 388 (2021).
  * Eckhoff  _et al._ (2020) M. Eckhoff, K. N. Lausch, P. E. Blöchl, and J. Behler, The Journal of Chemical Physics 153, 164107 (2020).
  * Gupta  _et al._ (2021) A. Gupta, S. Chakraborty, and R. Ramakrishnan, Machine Learning: Science and Technology 2, 035010 (2021).
  * von Rudorff  _et al._ (2020a) G. F. von Rudorff, S. Heinen, M. Bragato, and O. A. von Lilienfeld, Mach. Learn. Sci. Technol. 1, 45026 (2020a).
  * Bragato  _et al._ (2020) M. Bragato, G. F. von Rudorff, and O. A. von Lilienfeld, Chemical science 11, 11859 (2020).
  * Eckhoff and Behler (2021) M. Eckhoff and J. Behler, npj Computational Materials 7, 170 (2021).
  * Ko  _et al._ (2021) T. W. Ko, J. A. Finkler, S. Goedecker, and J. Behler, Accounts of Chemical Research 54, 808 (2021).
  * Huang and von Lilienfeld (2020b) B. Huang and O. A. von Lilienfeld, Nature Chem. 12, 945 (2020b).
  * Huang and von Lilienfeld (2020c) B. Huang and O. A. von Lilienfeld, “Dictionary of 140k gdb and zinc derived amons,” (2020c), arXiv:2008.05260 [physics.chem-ph] .
  * Goodall and Lee (2020) R. E. Goodall and A. A. Lee, Nature communications 11, 1 (2020).
  * Goodall  _et al._ (2021) R. E. A. Goodall, A. S. Parackal, F. A. Faber, R. Armiento, and A. A. Lee, “Rapid discovery of novel materials by coordinate-free coarse graining,” (2021), arXiv:2106.11132 [cond-mat.mtrl-sci] .
  * Wieder  _et al._ (2020) O. Wieder, S. Kohlbacher, M. Kuenemann, A. Garon, P. Ducrot, T. Seidel, and T. Langer, Drug Discovery Today: Technologies (2020), https://doi.org/10.1016/j.ddtec.2020.11.009.
  * Rupp  _et al._ (2012b) M. Rupp, A. Tkatchenko, K.-R. Müller, and O. A. von Lilienfeld, Phys. Rev. Lett. 108, 058301 (2012b).
  * Rogers and Hahn (2010) D. Rogers and M. Hahn, J. Chem. Inf. Mod. 50, 742 (2010).
  * Schütt  _et al._ (2018) K. T. Schütt, H. E. Sauceda, P.-J. Kindermans, A. Tkatchenko, and K.-R. Müller, The Journal of Chemical Physics 148, 241722 (2018).
  * Alibakhshi and Hartke (2021) A. Alibakhshi and B. Hartke, chemrxiv (2021), 10.26434/chemrxiv.14043503.
  * Riniker (2017) S. Riniker, Journal of chemical information and modeling 57, 726 (2017).
  * Axelrod and Gomez-Bombarelli (2020) S. Axelrod and R. Gomez-Bombarelli, “Molecular machine learning with conformer ensembles,” (2020), arXiv:2012.08452 [cs.LG] .
  * Gebhardt  _et al._ (2020) J. Gebhardt, M. Kiesel, S. Riniker, and N. Hansen, J. Chem. Inf. Mod. 60, 5319 (2020).
  * Scheen  _et al._ (2020) J. Scheen, W. Wu, A. S. J. S. Mey, P. Tosco, M. Mackey, and J. Michel, J. Chem. Inf. Mod. 60, 5331 (2020).
  * Lim and Jung (2021) H. Lim and Y. Jung, Journal of Cheminformatics 13 (2021), 10.1186/s13321-021-00533-z.
  * Vermeire and Green (2021) F. H. Vermeire and W. H. Green, Chemical Engineering Journal 418, 129307 (2021).
  * Venkatraman  _et al._ (2018) V. Venkatraman, S. Evjen, H. K. Knuutila, A. Fiksdahl, and B. K. Alsberg, Journal of Molecular Liquids 264, 318 (2018).
  * Kamel  _et al._ (2018) M. Kamel, C. Grulke, R. Judson, and A. Williams, Journal of Cheminformatics 10 (2018), 10.1186/s13321-018-0263-1.
  * Zaverkin  _et al._ (2022) V. Zaverkin, J. Netz, F. Zills, A. Köhn, and J. Kästner, Journal of Chemical Theory and Computation 18, 1 (2022).
  * Reinhardt and Cheng (2021) A. Reinhardt and B. Cheng, Nature communications 12, 1 (2021).
  * Cheng  _et al._ (2021) B. Cheng, M. Bethkenhagen, C. J. Pickard, and S. Hamel, arXiv preprint arXiv:2103.09035 (2021).
  * Monserrat  _et al._ (2020) B. Monserrat, J. G. Brandenburg, E. A. Engel, and B. Cheng, arXiv preprint arXiv:2006.13316 (2020).
  * Vapnik (1998) V. N. Vapnik,  _Statistical Learning Theory_ (Wiley-Interscience, 1998).
  * Deisenroth  _et al._ (2020) M. P. Deisenroth, A. A. Faisal, and C. S. Ong,  _Mathematics for Machine Learning_ (Cambridge University Press, 2020).
  * Schütt  _et al._ (2020) K. T. Schütt, S. Chmiela, O. A. von Lilienfeld, A. Tkatchenko, K. Tsuda, and K.-R. Müller,  _Machine Learning Meets Quantum Physics_ (Springer, 2020).
  * Koliński  _et al._ (2004) A. Koliński _et al._ , Acta Biochimica Polonica 51 (2004).
  * Liwo  _et al._ (2014) A. Liwo, M. Baranowski, C. Czaplewski, E. Gołaś, Y. He, D. Jagieła, P. Krupa, M. Maciejczyk, M. Makowski, M. A. Mozolewska, _et al._ , Journal of molecular modeling 20, 1 (2014).
  * Levitt and Warshel (1975) M. Levitt and A. Warshel, Nature 253, 694 (1975).
  * Becke (1988) A. D. Becke, Phys. Rev. A 38, 3098 (1988).
  * Perdew (1986) J. P. Perdew, Phys. Rev. B 33, 8822 (1986).
  * Duarte Ramos Matos  _et al._ (2017) G. Duarte Ramos Matos, D. Y. Kyu, H. H. Loeffler, J. D. Chodera, M. R. Shirts, and D. L. Mobley, J. Chem. Engin. Data 62, 1559 (2017).
  * Gibbs (1873) J. W. Gibbs, Transactions of the Connecticut Academy of Arts and Sciences 2, 382–404 (1873).
  * Fogolari  _et al._ (2002) F. Fogolari, A. Brigo, and H. Molinari, J. Mol. Recogn. 15, 377 (2002).
  * Mennucci  _et al._ (2002) B. Mennucci, J. Tomasi, R. Cammi, J. R. Cheeseman, M. J. Frisch, F. J. Devlin, and P. J. Stephens, J. Phys. Chem. A 106, 6102 (2002).
  * Marenich  _et al._ (2013) A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J. Chem. Th. Comp. 9, 609 (2013).
  * Marenich  _et al._ (2009) A. V. Marenich, C. J. Cramer, and D. G. Truhlar, J. Phys. Chem. B 113, 6378 (2009).
  * Klamt and Schüürmann (1993) A. Klamt and G. Schüürmann, J. Chem. Soc., Perkin Trans. 2 , 799 (1993).
  * Klamt (1995a) A. Klamt, J. Phys. Chem. 99, 2224 (1995a).
  * Klamt and Eckert (2000) A. Klamt and F. Eckert, Fluid Phase Equilibria 172, 43 (2000).
  * Beglov and Roux (1997) D. Beglov and B. Roux, J. Phys. Chem. B 101, 7821 (1997).
  * Kovalenko and Hirata (1998) A. Kovalenko and F. Hirata, Chem. Phys. Lett. 290, 237 (1998).
  * Roy and Kovalenko (2019) D. Roy and A. Kovalenko, J. Phys. Chem. A 123, 4087 (2019).
  * Rappoport and Furche (2010) D. Rappoport and F. Furche, The Journal of Chemical Physics 133, 134105 (2010).
  * Schäfer  _et al._ (1994) A. Schäfer, C. Huber, and R. Ahlrichs, The Journal of Chemical Physics 100, 5829 (1994).
  * Weigend and Ahlrichs (2005) F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005).
  * Weigend (2006) F. Weigend, Phys. Chem. Chem. Phys. 8, 1057 (2006).
  * Müller  _et al._ (1996) K.-R. Müller, M. Finke, N. Murata, K. Schulten, and S.-i. Amari, Neural computation 8, 1085 (1996).
  * Cortes  _et al._ (1993) C. Cortes, L. D. Jackel, S. A. Solla, V. N. Vapnik, and J. S. Denker, in  _NIPS_ (1993).
  * Cortes  _et al._ (1994) C. Cortes, L. D. Jackel, S. A. Solla, V. Vapnik, and J. S. Denker, in  _Advances in Neural Information Processing Systems_ (1994) pp. 327–334.
  * Kromann (2021) J. C. Kromann, “Calculate root-mean-square deviation (rmsd) of two molecules using rotation,”  (2021).
  * Walker  _et al._ (1991) M. W. Walker, L. Shao, and R. A. Volz, CVGIP: Image Understanding 54, 358 (1991).
  * Kabsch (1976) W. Kabsch, Acta Crystallographica Section A 32, 922 (1976).
  * Ruddigkeit  _et al._ (2012) L. Ruddigkeit, R. van Deursen, L. C. Blum, and J.-L. Reymond, Journal of Chemical Information and Modeling 52, 2864 (2012).
  * Zhao and Truhlar (2008) Y. Zhao and D. Truhlar, Theoretical Chemistry Accounts 120, 215 (2008).
  * Frisch  _et al._ (2016) M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman, and D. J. Fox, “Gaussian16 Revision C.01,” (2016), gaussian Inc. Wallingford CT.
  * Klamt (1995b) A. Klamt, The Journal of Physical Chemistry 99, 2224 (1995b).
  * Klamt  _et al._ (1998) A. Klamt, V. Jonas, T. Bürger, and J. C. Lohrenz, The Journal of Physical Chemistry A 102, 5074 (1998).
  * Eckert and Klamt (2002) F. Eckert and A. Klamt, AIChE Journal 48, 369 (2002).
  * Chung  _et al._ (2021) Y. Chung, F. H. Vermeire, H. Wu, P. J. Walker, M. H. Abraham, and W. H. Green, ChemRxiv (2021), 10.33774/chemrxiv-2021-djd3d.
  * Lipparini  _et al._ (2013) F. Lipparini, B. Stamm, E. Cancè, Y. Maday, and B. Mennucci, Journal of Chemical Theory and Computation 9, 3637–3648 (2013).
  * Sun  _et al._ (2018) Q. Sun, T. C. Berkelbach, N. S. Blunt, G. H. Booth, S. Guo, Z. Li, J. Liu, J. D. McClain, E. R. Sayfutyarova, S. Sharma, S. Wouters, and G. K.-L. Chan, WIREs Computational Molecular Science 8, e1340 (2018).
  * Perdew  _et al._ (1996) J. P. Perdew, M. Ernzerhof, and K. Burke, Journal of Chemical Physics 105, 9982 (1996).
  * Adamo and Barone (1999) C. Adamo and V. Barone, The Journal of Chemical Physics 110, 6158 (1999).
  * Onufriev  _et al._ (2004) A. Onufriev, D. Bashford, and D. A. Case, Proteins 55, 383 (2004).
  * Weiser  _et al._ (1999) J. Weiser, P. S. Shenkin, and W. C. Still, J. Comp. Chem. 20, 217 (1999).
  * Roy and Kovalenko (2021) D. Roy and A. Kovalenko, J 4, 604 (2021).
  * Kovalenko and Hirata (1999) A. Kovalenko and F. Hirata, The Journal of Chemical Physics 110, 10095 (1999).
  * Reimann and Kaupp (2020) M. Reimann and M. Kaupp, The Journal of Physical Chemistry A 124, 7439 (2020).
  * Chung  _et al._ (2020) Y. L. Chung, R. J. Gillis, and W. H. Green, Aiche Journal 66 (2020).
  * Chung  _et al._ (2022) Y. Chung, F. H. Vermeire, H. Wu, P. J. Walker, M. H. Abraham, and W. H. Green, Journal of Chemical Information and Modeling 62, 433 (2022).
  * Zhang  _et al._ (2022) D. Zhang, S. Xia, and Y. Zhang, Journal of Chemical Information and Modeling 62, 1840 (2022).
  * Ramakrishnan  _et al._ (2015) R. Ramakrishnan, P. O. Dral, M. Rupp, and O. A. von Lilienfeld, J. Chem. Th. Comp. 11, 2087 (2015).
  * Nigam  _et al._ (2021) A. Nigam, R. Pollice, M. Krenn, G. d. P. Gomes, and A. Aspuru-Guzik, Chem. Sci. 12, 7079 (2021).
  * Wu  _et al._ (2018) Z. Wu, B. Ramsundar, E. Feinberg, J. Gomes, C. Geniesse, A. S. Pappu, K. Leswing, and V. Pande, Chem. Sci. 9, 513 (2018).
  * Gao  _et al._ (2021) P. Gao, X. Yang, Y.-H. Tang, M. Zheng, A. Andersen, V. Murugesan, A. Hollas, and W. Wang, Phys. Chem. Chem. Phys. 23, 24892 (2021).
  * Faber  _et al._ (2018) F. A. Faber, A. S. Christensen, B. Huang, and O. A. von Lilienfeld, J. Chem. Phys. 148, 241717 (2018).
  * Christensen  _et al._ (2020) A. S. Christensen, L. A. Bratholm, F. A. Faber, and O. A. von Lilienfeld, J. Chem. Phys. 152, 044107 (2020).
  * Moussa (2012) J. E. Moussa, Phys. Rev. Lett. 109, 059801 (2012).
  * Huang and von Lilienfeld (2016b) B. Huang and O. A. von Lilienfeld, The Journal of Chemical Physics 145, 161102 (2016b).
  * Parsaeifard  _et al._ (2020) B. Parsaeifard, D. De, A. Christensen, F. Faber, E. Kocer, S. De, J. Behler, A. Lilienfeld, and S. Goedecker, Machine Learning: Science and Technology 2 (2020), 10.1088/2632-2153/abb212.
  * Langer  _et al._ (2022) M. Langer, A. Goeßmann, and M. Rupp, npj Computational Materials 8, 41 (2022).
  * Pozdnyakov  _et al._ (2020) S. N. Pozdnyakov, M. J. Willatt, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, Phys. Rev. Lett. 125, 166001 (2020).
  * von Lilienfeld  _et al._ (2015) O. A. von Lilienfeld, R. Ramakrishnan, M. Rupp, and A. Knoll, International Journal of Quantum Chemistry 115, 1084 (2015), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.24912 .
  * Heinen  _et al._ (2021) S. Heinen, G. F. von Rudorff, and O. A. von Lilienfeld, The Journal of Chemical Physics 155, 064105 (2021).
  * von Rudorff  _et al._ (2020b) G. F. von Rudorff, S. Heinen, M. Bragato, and O. A. von Lilienfeld, Mach. Learn. Sci. Technol. 1, 45026 (2020b).
  * Parsaeifard  _et al._ (2021b) B. Parsaeifard, D. S. De, A. S. Christensen, F. A. Faber, E. Kocer, S. De, J. Behler, O. A. von Lilienfeld, and S. Goedecker, Machine Learning: Science and Technology 2, 015018 (2021b).
  * Behler (2011) J. Behler, The Journal of chemical physics 134, 074106 (2011).
  * Bircher  _et al._ (2021) M. P. Bircher, A. Singraber, and C. Dellago, Machine Learning: Science and Technology 2, 035026 (2021).
  * Pracht  _et al._ (2020) P. Pracht, F. Bohle, and S. Grimme, Phys. Chem. Chem. Phys. 22, 7169 (2020).
  * Grimme (2019) S. Grimme, Journal of Chemical Theory and Computation 15, 2847 (2019).
  * Pracht and Grimme (2021) P. Pracht and S. Grimme, Chem. Sci. 12, 6551 (2021).
  * Yoshikawa and Hutchison (2019) N. Yoshikawa and G. Hutchison, Journal of Cheminformatics 11 (2019).
  * Riniker and Landrum (2015) S. Riniker and G. A. Landrum, Journal of Chemical Information and Modeling 55, 2562 (2015).
  * Miteva  _et al._ (2010) M. Miteva, F. Guyon, and P. Tufféry, Nucleic Acids Research 38, W622 (2010).
  * Hawkins  _et al._ (2010) P. C. D. Hawkins, A. G. Skillman, G. L. Warren, B. A. Ellingson, and M. T. Stahl, Journal of Chemical Information and Modeling 50, 572 (2010).
  * Landrum  _et al._ (2006) G. Landrum _et al._ , software (2006).
  * Boyle  _et al._ (2011) N. M. O. Boyle, T. Vandermeersch, C. J. Flynn, A. R. Maguire, and G. R. Hutchison, Journal of Chemical Information and Modeling 3 (2011), 10.1186/1758-2946-3-8.
  * Vainio and Johnson (2007) M. J. Vainio and M. S. Johnson, Journal of Chemical Information and Modeling 47, 2462 (2007).
  * Mansimov  _et al._ (2019) E. Mansimov, O. Mahmood, S. Kang, and K. Cho, Scientific Reports 9 (2019).
  * Kästner  _et al._ (2009) J. Kästner, J. Carr, T. W. Keal, W. Thiel, A. Wander, and P. Sherwood, The journal of physical chemistry. A 113 43, 11856 (2009).
  * Chan  _et al._ (2018) L. Chan, G. Hutchison, and G. Morris, chemrxiv (2018), 10.26434/chemrxiv.7228940.
  * Meldgaard  _et al._ (2021) S. A. Meldgaard, J. Köhler, H. L. Mortensen, M.-P. V. Christiansen, F. Noé, and B. Hammer, “Generating stable molecules using imitation and reinforcement learning,” (2021), arXiv:2107.05007 [physics.chem-ph] .
  * Wu  _et al._ (2020) H. Wu, J. Köhler, and F. é, arxiv:2002.06707 (2020).
  * Jakalian  _et al._ (2000) A. Jakalian, B. L. Bush, D. B. Jack, and C. I. Bayly, J. Comp. Chem. 21, 132 (2000).
  * (163) “TURBOMOLE V7.2 2017, a development of University of Karlsruhe and Forschungszentrum Karlsruhe GmbH, 1989-2007, TURBOMOLE GmbH, since 2007; available from
http://www.turbomole.com.” .

  * (164) “Cosmoconf 4.3,” COSMOlogic GmbH Co. KG; http://www.cosmologic.de, Leverkusen, Germany.
  * Eckert and Klamt (2018) F. Eckert and A. Klamt, “Cosmotherm,” (2018), bIOVIA COSMOtherm, Release 2021; Dassault Systèmes. http://www.3ds.com.
  * Lemm  _et al._ (2021b) D. Lemm, G. F. von Rudorff, and A. von Lilienfeld, “Leruli.com, online molecular property predictions in real time and for free,” www.leruli.com (2021b).